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ABSTRACT 

A  complete  description  is  given  of  the  design,  implementation  and  use  of  a 
family  of  very  fast  and  efficient  large  scale  minimum  cost  primal  network  programs. 
Choice  of  data  structures  and  computational  testing  of  the  network  system  GNET  are 
discussed.   Important  extensions  are  explained  such  as  exploitation  of  special 
problem  structure,  element  generation  techniques,  post  optimality  analysis,  operation 
with  problem  generators  and  external  problem  files,  and  generalization  beyond  pure 
network  models. 


INTRODUCTION 

This  paper  reports  the  development  of  a  large  scale  primal  network  code  that  is 
possibly  the  fastest  and  most  efficient  program  currently  available  for  solving 
capacitated  transshipment  problems.   The  capacitated  transshipment  problem  is  the 
most  general  of  the  minimum  cost  flow  models  which  include  the  capacitated  and 
uncapacitated  transportation  problems  and  the  personnel  assignment  problem.   These 
models  are  used  for  a  large  number  of  diverse  applications  that  include  transporta- 
tion of  goods,  design  of  communications  and  pipeline  systems,  assignment  of  men  to 
jobs,  bid  evaluation  and  production  planning.   For  further  discussion  of  these 
applications  see  the  survey  articles  of  Bradley  [5],  Elmaghraby  [17]  and  Fulkerson 
[21]  and  the  textbooks  of  Busacher  and  Saaty  [8],  Charnes  and  Cooper  [9],  Dantzig 
[14]  and  Ford  and  Fulkerson  [19]. 

The  capacitated  transshipment  model  and  its  specializations  are  minimum  cost 

network  flow  problems.   The  goal  is  to  determine  how  (or  at  what  rate)  a  good  should 

flow  through  the  arcs  of  a  network  to  minimize  shipment  costs.   The  network  is  a 

directed  graph  defined  by  a  set  of  nodes,  hi,      and  a  set  of  arcs,   A,   with  ordered 

pairs  of  nodes  (tail,  head)  as  elements  indexed  by  k.   For  each  arc  there  is  a 

shipping  cost  per  unit  flow,   c,  ,   a  minimum  allowable  flow  (or  lower  bound),  &,, 

and  a  maximum  allowable  flow  (or  capacity),  a,   .   Each  node  is  either  a  supply  node 

where  units  of  the  good  enter  the  network,  a  demand  node  where  units  leave,  or  a 

transshipment  node.   The  problem  is  to  minimize  total  costs  with  flows,   x,  ,   that 

satisfy  the  associated  lower  bounds  and  capacities  and  preserve  the  conservation  of 

flow  at  each  node:  min   J   ci,xi, 

keA 

s.t.  I  xk     -  I  xk  =  bif  i   e   W  (1) 

keA  with  keA  with 

tail  i  head  i 


£^x,^6L  keA 


where  b.  =  {  supply  if  i  is  a  supply  node;  -demand  if   i  is  a  demand  node;  0 
otherwise  }. 

We  choose  this  notation  to  emphasize  that  data  will  be  stored  only  for  arcs  that 

are  present  in  the  network,  and  operations  are  defined  only  for  use  with  these  arcs. 


The  usual  textbook  notation  with,  for  instance,   £  I   x    is  particularly  inappropriat 

i  J 
since  for  large  practical  problems  it  is  never  true  that  all  node  pairs  are  connected 

by  an  arc.   Further,  our  notation  allows  multiple  arcs  (connecting  any  given  pair  of 

nodes)  which  are  common  in  practical  problems.   This  notation  is  also  consistent  with 

the  input  requirements  of  all  contemporary  large  scale  network  optimization  and  lineal 

programming  codes. 

These  models  can  be  solved  as  linear  programming  problems  with  a  constraint  for 
each  node  and  a  variable  for  each  arc.   For  large  scale  problems,  contemporary 
commercial  linear  programming  codes  require  50-200  times  as  much  computer  time  and 
considerably  more  space  for  data  storage  than  special  purpose  network  flow  algorithms. 

The  transportation  model  was  originally  proposed  by  Hitchcock  [36]  1941  and 
Koopmans  {40]  1946.   Both  presented  computational  methods  that  would  now  be  called 
"primal  simplex."  Hitchcock  showed  that  an  optimal  solution  would  be  an  extreme 
point  solution  and  showed  how  to  iteratively  construct  better  extreme  point  solutions. 
He  noted  alternate  optimal  solutions  and  degeneracy  (all  in  6  1/4  pages) ,   Koopmans 
developed  simplex  multipliers,  or  "node  potentials"  and  the  optimality  criterion  and 
he  showed  that  an  extreme  point  is  equivalent  to  a  tree.   Dantzig  [13]  1951  showed 
how  the  transportation  problem  could  be  solved  by  his  simplex  algorithm;  he  also 
developed  a  special  variant  of  the  simplex  algorithm  for  the  transportation  problem. 
Orden  [48]  showed  that  these  results  can  be  extended  to  the  transshipment  problem. 

Network  models  are  widely  used  because  they  accurately  describe  a  variety  of 
important  applications,  and  network  algorithms  can  efficiently  solve  large  problems 
(many  thousands  of  equations  and  variables) .   Such  models  are  popular  because  they 
are  more  readily  accepted  by  nonanalysts  than  perhaps  any  other  class  of  opera- 


tions  research  models.   Network  algorithms  can  economically  solve  problems  with  more 
variables  than  virtually  any  other  optimization  technique.   There  has  been  a  surge  of 
interest  in  these  models  because  more  efficient  computer  programs  have  made  possible 
the  economic  solution  of  much  larger  problems. 

Although  many  papers  have  been  written  in  this  general  area,  and  significant 
computational  breakthroughs  have  been  reported,  there  has  not  previously  been  a  single, 
unified  description  of  a  complete  implementation,  nor  have  "new  generation"  computer 
programs  been  made  generally  available.   Here  we  report  the  research  and  computational 
experiments  which  have  produced  GNET,  an  extremely  efficient  yet  relatively  simple 
code.   An  important  objective  of  this  paper  is  to  make  these  new  approaches  easily 
accessible  to  a  wide  audience  via  a  clear  mathematical  exposition  and  a  concrete 
Example  of  a  highly  efficient  FORTRAN  program.   Further,  the  availability  of  the 
computer  program  will  now  make  it  possible  for  any  investigator  to  reproduce  and 
extend  the  experimental  results. 

The  approach  we  have  used  in  this  research  is  to  study  the  data  structures  which 
seem  to  be  most  fundamental  in  the  sense  that  they  can  be  applied  to  many  types  of 
mathematical  programming  situations.   In  this  context,  we  view  the  major  advances 
over  the  last  thirty  years  in  efficient  solution  of  large  linear  programming  problems 
(for  example:   revised  simplex,  product  form  inverse,  bounded  variables,  generalized 
upper  bounding,  factorization,  sparse  matrix  methods,  etc.)  as  major  changes  in  the 

representation  of  the  data  accompanied  only  by  necessary  modifications  of  the  simplex 

I 

procedure  to  accommodate  these  new  data  structures.   The  computational  breakthroughs 
in  primal  network  codes  are  also  due  to  improvements  in  data  representation  and  renew 
our  interest  in  the  subtle  relationships  between  the  algorithms  and  the  data  struc- 
tures employed  for  implementation. 

It  is  helpful  to  view  networks  as  an  important  special  case  of  large  scale  linear 
programming.   This  approach  is  crucial  and  overdue  because  there  has  been  very  little 


success  in  extending  the  graph  theoretic  basis  tree  approach  to  more  general  mathe- 
matical programming  models.   For  instance,  the  graph  theoretic  proofs  of  the 
mathematical  correctness  of  generalizations  of  pure  network  minimization  problems  are 
arduous.   Thus,  in  this  respect  we  consider  a  purely  graph  theoretic  approach 
to  be  a  dead  end.  We  have  taken  an  entirely  different  tack  in  developing  the  mathe- 
matical framework  of  a  general  linear  programming  problem  and  specializing  it  to 
primal  network  models.   Therefore,  it  is  necessary  to  build  a  sufficient  mathematical 
foundation  to  answer  the  question:   "What  is  the  capacitated  transshipment  problem  a 
specialization  of?",  rather  than,  "How  can  we  generalize  basis  trees?"  Thus  we  invit 
the  reader  to  view  this  paper  as  drawing  from  a  general  theory  of  large  scale  mathe- 
matical programming  for  which  the  network  algorithms  are  concrete  realizations  of  a 
rich  algebraic  view  of  a  problem  with  special  structure. 

Approaches  other  than  the  primal  simplex  have  been  proposed  including  out-of- 
kilter  (Fulkerson  [20]),  primal-dual  (Ford  and  Fulkerson  [18]),  dual  (Balas  and  Hamme 
[1])  ,  path  (Busacker  and  Gowen  [7]),  negative  cycle  (Klein  [38])  and  scaling  (Edmonds 
and  Karp  [16]).   Also,  special  algorithms  have  been  developed  for  the  assignment, 
shortest  route  and  maximum  flow  problems. 

The  contemporary  implementations  of  the  primal  simplex  network  algorithm  are 
based  on  compact  representation  of  the  basis,  determination  of  the  outgoing  arc  withe 
trial  and  error  and  efficient  techniques  to  update  the  simplex  multipliers  at  each 
pivot.   Some  of  these  developments  were  suggested  in  the  1960 's.   Glicksman,  L. 
Johnson  and  Eselson  [25]  described  a  transportation  problem  with  few  sources  and  man)) 
sinks;  their  data  structure  may  be  viewed  as  storing  the  basis  as  an  arborescence  anc 
using  this  structure  to  efficiently  find  the  outgoing  arc.   E.  Johnson  [37]  describee 
a  triple  label  scheme  that  represents  the  basis  as  an  arborescence  and  allows  the 
simplex  multipliers  to  be  efficiently  updated.   Johnson  describes  his  work  as  a  modi- 
fication of  the  proposal  of  Scoins  [51]. 


The  proposals  of  Glicksman,  L.  Johnson  and  Eselson  [25],  Scoins  [51]  and 
E.  Johnson  [37]  for  primal  algorithms  were  not  immediately  pursued;  the  most  widely 
known  code  of  the  decade  was  an  out-of -kilter  implementation  by  Clasen  [12],   The 
contemporary  work  on  network  algorithms  was  begun  in  1970  by  Srinivasan  and  Thompson 

[54,  55],  Glover,  Karney  and  Klingman  [26]  and  Glover,  Karney,  Klingman  and  Napier 

[27].   This  work  was  a  break  with  the  past  in  that: 

1.  Primal  algorithms  were   considered  despite   all  the   experiments   in   the    1950' s   and 
early   1960's   that   showed   the   apparent   superiority   of   the   out-of-kilter  algorithm, 

2.  Contemporary   computer  science   tools   that  had  not  been   available   a  decade   earlier 
were   used,    and 

3.  Computer   codes  were   developed   for  much   larger  problems. 

Later  and  independently,   McBride    [44]    and  Graves   and  McBride    [32]    specialized   their 
tfork   on   factorization   of    linear  programs   to   transshipment   problems.      Although   their 
development  was   quite   different,    the  network  specialization   of   their   data  structures 
is   similar  in  many   respects    to  data  structures   that   evolved   from  a  graph  theoretic 
view  of  networks.      Mulvey    [45]   has    developed  an  efficient    large   scale  primal   code   at 
TRW.      Harris    [33]   has    developed   a  primal   algorithm   for   transportation  problems  with 
many   sinks    and   few  sources.      Langley,   Kennington   and  Shetty    [42]   have   also   developed 
a  primal   transshipment    code.. 

A  significant   aspect    of   contemporary  network   research  has  been   the   computational 
testing  of   different   algorithms   on   large   standard  test  problems.      One  major  topic  has 
been  primal  algorithms   versus    out-of-kilter  algorithms.      Experiments   in   the   1950' s   and 
early   1960's    convinced   researchers   that    the   out-of-kilter  algorithm  was   superior, 
especially   for  transshipment   problems.      The  most    comprehensive   recent    comparison  has 
been   done  by  Glover,   Karney   and  Klingman    [26]    and  Barr,    Glover  and  Klingman    [2]   who 
compare   the   algorithms   on   a  diverse   set   of   test   problems    [39].      Their  primal   code  was 
from  30%    (for  transshipment)    to   40%    (for  transportation)    faster  than   their  out-of-kilter 


code.   The  success  of  the  primal  algorithm  has  been  independently  verified  by  the 
experiments  of  others.   Most  researchers  now  believe  that  the  primal  algorithm  is 
superior  to  others  including  the  out-of-kilter  algorithm  (an  exception  is  Hatch  [35]) 
Current  primal  implementations  are  faster,  require  less  storage,  are  more  suitable  wh 
using  secondary  storage  devices  and  are  compatible  as  embedded  parts  of  more  general 
optimization  systems. 

Although  the  algebraic  development  of  our  computer  code  preceded  the  choice  of 
data  structures  and  the  algorithm,  we  can  also  establish  the  graph  theoretic  interpre 
tation  for  the  pure  network  problem.   Thus,  although  our  derivation  was  dissimilar 
from  historical  work  in  this  field,  we  are  able  to  show  the  relationship  of  our  work 
to  that  of  others.   For  expository  purposes  we  will  draw  from  both  linear  algebra  and 
graph  theory,  using  pictures  and  terminology  consistent  with  past  literature;  the 
detailed  algebraic  development  is  given  in  a  companion  report. 

We  continue  now  with  a  brief  algebraic  description  of  the  general  bounded 
variable  simplex  algorithm  and  several  commonly  used  implementation  options.   The 
algebraic  specialization  of  the  simplex  method  for  pure  network  problems  is  presented 
After  this  necessary  but  somewhat  mathematically  involved  section,  the  specific  desig 
decisions  and  experiments  carried  out  with  GNET  are  described,  including  computationa 
evidence  which  indicates  that  the  code  produced  by  this  approach  compares  very  favor- 
ably with  other  algorithms,  proprietary  or  otherwise.   Several  extensions  of  GNET  are 
presented,  including  codes  tailored  to  capacitated  and  uncapacitated  transportation 
problems,  and  other  variants  to  exploit  special  network  structure.   Postoptimal  and 
reoptimization  procedures  using  GNET  are  discussed.   Finally,  a  review  of  the  litera- 
ture traces  the  original  contributions  found  to  be  of  fundamental  importance  in  this 
work. 

THE  PRIMAL  SIMPLEX  ALGORITHM 
In  this  section  we  briefly  review  the  mathematical  underpinnings  of  the  bounded 
variable  revised  simplex  method.   In  order  to  maintain  a  broad  scope,  the  presentatio 
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intentionally  avoids  tangential  issues  of  specific  methods  and  implementations. 
Rather,  the  class  of  algebraic  algorithms  is  characterized  with  only  the  briefest 
indication  of  options  often  employed  for  actual  linear  programming  codes.   Small 
Roman  letters  will  denote  column  vectors,  and  a  prime  indicates  transpose.   Large 
letters  denote  matrices;  superscripts  denote  columns  and  subscripts  denote  rows. 

Consider 

min   c'x  s.t.   Ax  =  b   and  0  ^  x  ^  U.  ;  (2) 

where  A  is  viewed  for  the  present  as  a  matrix  of  technological  coefficients. 

As  a  practical  matter,  lower  bounds  on  the  variables  in  (1)  have  been  elimi- 
nated by  transformation. 

The  upper  bounds,  u,      are  most  readily  accommodated  implicitly.   Whenever  x 

k 

reaches  its  upper  bound,  it  is  logically  replaced  (fiO-fatO-Ctad)   by  u,  -  x,  ;   column  k 
is  logically  treated  as  if  its  sign  has  changed  and  the  explicit  right  hand  side  is 
transformed  to  b  -  A  U.  .   If  a  record  is  kept  of  each  variable  that  is  at  its  upper 
bound,  the  original  problem  solution  is  easily  recovered. 

By  construction  (possibly  including  introduction  of  unit  vectors  representing 
slack  and  artificial  variables)   A  may  be  partitioned  A  =  (B,N),  with  B   an  m  x  m 
matrix  of  linearly  independent  columns  which  is  called  a  basis. 

Given  a  feasible  basis  there  always  exists  a  unique   x  such  that 

Bx  =  b  .  (3) 


In  terms  of  this   x   there  is  always  a  current  basic  solution  x°  = 
partitioning   c   in  the  same  manner  as   A,   one  has 


and,    upon 


C'X°    :      (CB'CN}     10 


=    CBX 


(4) 


Any   generic  solution      x     satisfying  the   constraints    can  be   rewritten   in   terms   of   the 
basic  solution: 


Ax  =    (B,N) 


V 


=  BxB  +  NXn  =  b    . 


(5) 


Further,  since  B  is  a  basis,  there  exists  a  unique  transformation  Z  such  that 

BZ  =  N  .  (6) 

Multiplying  (6)  by  x^  and  subtracting  from  (3)  yields 

B(x-ZxN)  +  N^  =  b  ,  (7) 

and  subtracting  (7)  from  (5)  yields 

BCXg-lS-Z^J)  =  0   . 

Since  the  columns  of  B   are  linearly  independent,   x  =  x  -  Zx^   and  the  general 
solution  becomes 


x  = 


V 


(x_zv 


(8) 


With  this  form  it  is  easy  to  compare  the  value  of  x  to  any  current  solution 
x°   and  identify  an  improved  solution  when  one  exists.   The  value  of  the  generic 
solution  (use  (4)  and  (6))  is 

C'X  =  ci*B  +  ^  =  CB(*-ZXN}  +  ^ 

"  CB*  +  (CN"CBZ)XN 

=  c^x  +  (CjJ-u'N)^  ;  (9) 

where  u   (often  called  the  dual  solution  or  simplex  multipliers) is  the  solution  of 

u'B  =  cl    .  (10 

D 

From  (  9 )  it  is  clear  (since  x^  ^  0)  that  a  necessary  condition  for  an  improved 
solution  is  that  there  exist  a  column  of  N,   N  ,   such  that 

c^  -  u'Nk  <  0  .  (11 

With  such  a  column  consider  a  specific  vector  obtained  from  (8)  by  taking  all 
components  of   x^   zero  except   x,  , 


x'  =  (-x_z  ^V0'** ' ,xk,•  •  • »°) 


(12 


As  a  function  of  x^  ^  0,   this  vector  must  by  its  derivation  satisfy  the 
explicit  constraints,   Ax  =  b ;  feasibility  also  requires  satisfaction  of  the  bounds, 
0  ^  x  ^  Li.      For  components   z   >  0   this  requires   x.  -  z.i*.  ^  0,   or 

\  *  Vzik  •  <13) 

For  components   z   <  0   this  requires   x.  -  z .,  X  ^  ii.  ,      or 

^^  (x.-u.)/-Z.k  .  (14) 

For  V  \  s   V 

If   the    least  bound  on     x,       is      ll  ,      then     x,       stays   out   of   the  basis  but   goes   to 

its   upper  bound   and    (5)    and    (12)    yield     x  =   x  -  Z  it       as   a  new  basic  solution  with 
Bx  =  b   -  Nk6L  . 

If    (13)   is    the    least  bound   on     x,  ,      taking     x^   =   x./z  with      z        >   0      in    (12) 

i  k 

leads    to    the  exchange   of     B        and     N        in   the  partition  between  basic  and  nonbasic 

columns   and   the  new  basic  solution 

x     =   x     -   z   ,  (x.  /z..  )  ,  r  ^   i    ; 

r  r  rk     l      lk 

and     x.    =   x./z.,     . 
i  l      lk 

If    (14)    is   the    least  bound   on      x  ,      taking     x^   =    (u.-x.)/-z  with      z.      <   0     in 

i  k 

(12)    requires    the  basis   exchange   of      B        and     N        and  yields      x     =   x     -   z  ,  (u..-x. ) /-z ..  , 

rrrkiiik 

r  ^   i;      and     x.    =    (u.-x.)/-z..  ;      as   a  new  basis  with     x.      at  its   upper  bound, 
l  l      l  lk  l 

Assume    that    there   is    a  current  basis      B,      a  current   solution     x     to     Bx  =  b ,      and 
a   current   solution      u     of      u'B  =   c' .      A  step   of   the   simplex  procedure   is   summarized: 

D 

k 

51.  VnLoJLOVJL.      Select   a   candidate   variable   to   enter  the  basis  with      (c, -u'N   )    <   0. 

k  k 

52.  Ratto   Tut.      Find  the   greatest  bound  such   that    (with     BZ     =  N   ): 

a)  ^  *  V 

b)  \*  VZrk      f°r      Zrk   *   °» 

c)  x,    ^    (a  -x   )/-z   .       for      z   .     <   0. 
It  r     r  rk  rk 


If  the  minimum  ratio  is  determined  by  case  b)  or  c)  ,  let   i  be  a  basis  variable  for 
which  the  minimum  is  achieved. 

S3.  ?AJV0t.      Update  the  primal  solution:   x  =  x  -  x,N  .   If  x,   is  bounded  by  case 

a),  reflect  x   and  leave  the  basis  and  dual  solution  unchanged.   Otherwise, 
k 

i        k 
change  the  basis  by  exchanging  B   and  N  ,  for  case  c)  reflect   x. .   For  case 

b)  or  c)  find  the  new  dual  solution  to  u'B  =  c '  . 

In  executing  the  simplex  algorithm  a  number  of  options  have  customarily  been 

k    k 
employed  for  generating  the  solutions  of  the  linear  systems   BZ  =  N   and  u'B  =  c'. 

D 

In  general  algorithms,  the  basis  inverse  Q  =  B    is  usually  used,  stored  and 
updated  in  some  form.   Further,  although  there  is  no  difficulty  in  deriving  a  new 
algebraic  solution  to  (10),  u'B  =  c'  ,   as  a  practical  matter  u  may  be  directly 
achieved  from  u  by  simple  transformation. 

Proposition:   u  =  u  +  AQ.   where   Q.   is  the  i —  row  of  the  inverse  of  B.       (II 

k  k     k        k      k         I 

The  new  basic  column  N    determines   A   as   uN  =  uN  +  A(Q.N  )  =  uN  +  A(Q.BZ 

k  k  i 

uN  +  ^z-\<   ~   %    '   so  tnat   ^  =  (c,-uN  )/z.,  .   Exclusive  of  the  outgoing  column  B  , 

u'Br  =  uBr  +  A(Q.B  )  =  uB  =  c  ,      r^i.   n 
1  r      r    r 

i        k 
The  (pivotal)  update  of  Q  after  exchange  of  B   and  N   is  easily  derived 

(e.g.,  [30]).   The  most  elementary  and  explicit  procedure  is  to  carry  and  update  (by 

pivoting)    a   complete   tableau 

ON  Q> 

,vu'N      . 

The   revised  simplex  procedure   generates   these   elements   as   needed  by   access   to   columns 
of     N,      c     and     b      and  use   of     Q.      Most    full  scale   systems   employ   an  additional 
refinement  by  expressing     Q     as   the   product   of  elementary   "eta"   column  vectors,   each 
representing   the   pivotal   transformations   generating     Q      from  an  initial  basis. 
Frequently   the  history   of    'eta'    columns   grows   too   long  for   reasonably  efficient 
generation   of     Q,      or  numerical  error  is   propagated  and  detected,    forcing   a  rein- 
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version  with  'eta'  column  representatives  from  only  the  current  basis. 

Other  systems  support  the  solution  procedure  by  using  combinations  of  features 
such  as  an  "LU"  decomposition  of  B   {3],  a  Cholesky  decomposition  of  BB'  =  LL1 
[23,  49,  50],  or  list  representation  of  nonzero  elements  of  problem  components  and 
coefficient  representation  by  pointers  to  a  table  of  real  values.   Hybrid  schemes 
factor  B   into  partitions  with  special  structure:   Generalized  Upper  Bounding  (GUB) 
identifies  an  inherent  identity  matrix  for  some  rows  of  B   [15];  a  partial 
triangulation  of   B  with  an  inverse  for  remaining  columns  .and  rows  can  be  used  [31 
32].   Whether  systems  solve  (6)  and  (10)  by  triangular  substitution,  inverse  trans- 
formation, or  some  combination,  all  are  algebraically  equivalent  simplex  implemen- 
tations differing  only  in  the  structures  chosen  to  support  computation  for  the  class 
of  problems  at  hand. 

PRIMAL  NETWORK  SPECIALIZATION 

A  specialization  of  the  simplex  algorithm  to  the  transportation  problem  was 
developed  by  Dantzig  [13]  in  1951.   It  is  not  surprising  that  the  transportation 
algorithm  was  developed  immediately  after  the  simplex  algorithm,  because  the  works  of 
Hitchcock  [36]  1941  and  Koopmans  [40]  1946  on  the  transportation  problem  contain  many 
concepts  that  presage  the  simplex  algorithm.   The  interaction  of  general  linear 
programming  algorithms  and  transshipment  algorithms  has  a  long  history  that  has 
enriched  the  study  of  both.  . 

Here  we  establish  explicitly  the  relationship  between  the  general  primal  simplex 
algorithm  and  the  modern  implementations  of  the  transshipment  algorithms.   Our  goal  is 
to  understand  the  algebraic  foundations  of  the  modern  transshipment  implementations. 
Also  (and  perhaps  more  importantly)  we  lay  the  groundwork  for  the  next  stage  in  the 
interplay  of  these  models:   the  incorporation  into  the  next  generation  of  general 
linear  programming  computer  systems  of  the  important  ideas  that  have  made  possible  the 
breakthroughs  for  transshipment  problems. 

The  fundamental  fact  that  permits  design  of  efficient  primal  transshipment 
algorithms  is  the  well-known  result  that  any  transshipment  basis  can  be  put  in  (upper) 


triangular  form  by  a  simple  permutation  triangulation.   This  inherent  triangularity 
can  be  exploited  by  network  specializations  of  the  simplex  method  by  directly  solving 
(6)  by  back  substitution  and  (10)  by   forward  substitution.   Also,  the  triangulate1 
basis  simplex  algorithms  lead  to  much  more  efficient  network  solutions  due  to  other 
fortunate  simplifications.   The  most  remarkable  of  these  is  that  the  solution  update  \ 
of  step  S3  can  be  accompanied  (in  fact  aided)  by  a  very  simple  and  efficient  dynamic 
retriangulation  of  each  new  basis. 

An  initial  transshipment  basis  with  full  row  rank  can  always  be  constructed  by 
introducing  for  each  row  in  (1)  a  unit  vector  with  sign  matching  that  of  the  right 
hand  side  (negative  for  demand  nodes) .   With  the  addition  of  these  artificial  vectors 
A  has  full  row  rank  and  each  column  of  A  has  a  single  nonzero  element  1,  a  single 
nonzero  element  -1,  or  two  nonzero  elements  (a  1  and  a  -1). 

Theorem:   Any  basis  B  extracted  from  A  for  a  transshipment  problem  can  be 
triangulated  by  rearranging  rows  and  rearranging  columns. 

Proof:   Let  B  have  m  rows.   Locate  a  column  with  a  single  nonzero  entry.   Exchang 
rows  and  columns  so  that  the  nonzero  entry  is  the  first  diagonal  element.   For  the 
k —  step  of  the  construction  rows  and  columns  have  already  been  rearranged  so  that 
columns   l,2,...,k-l  have  only  zeros  below  the  diagonal.   Select  a  column  with  a 
single  nonzero  entry  in  rows   k,...,m.   Exchange  rows  and  columns  so  that  the  nonzero 
element  is  the  k —  diagonal  element.   This  construction  is  well  defined  only  if 
there  is  a  column  with  a  single  nonzero  element  in  rows  k,...,m.   There  must  be  such 
a  column  for  otherwise  each  column  would  have  no  nonzero  elements  or  exactly  one  +1 
and  one  -1  in  rows   k,...,m  and  the  sum  of  rows   k,...,m  would  be  the  zero  row  whic 
would  contradict  the  assumption  that  B   is  a  basis.    D 

A  graph  can  be  defined  that  represents  the  transshipment  basis.   Let   I  =  {i.  ,i„ 
•••>i  }  be  a  row  ordering  corresponding  to  a  triangulated  B.   Associate  with  each 
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lode      i,       the   row  number     p(i,  )      of  the  off  diagonal  element   in   the     k —     column  of   the 
k  k- 

liriangulated  basis;    if  there   is  no  nonzero   off  diagonal  element,    set     p(i,  )      to     m  +  1. 
)efine   a  graph  with  nodes      l,2,...,m+l     and  let      (i    ,p(i   ))      k  =   1,2,...  ,m     be   directed 
ires   from     i,       to     p(i   ).      Since   each  node      i        is    connected  to  a  node     p(i   )    =   i 
pith     h   <  k      or   to  node     m  +   1,    there   is    a  directed  path   from  each  node  to  node     m  +  1. 
'he   graph   is    called  a  HJOOt&d  CUibofieJ>C(inCL2.    [4]   with  node     m  +   1     the  KDOt.      Ignoring 
:he   orientation   of   the   arcs,    since   the   graph   is.  connected  and  has     m  +   1     nodes   and 
a     arcs   it   is    a   tree    (it    can  be   shown    [4]    that    this   definition   is   equivalent   to   the 
lsual  definition   of   a   tree   as   a   connected  graph  with  no   cycles) .      In   the   computer 
cience   literature    (e.g.,    [41])    the    term  tree   is   often   used  instead  of   rooted 
irborescence.      Figure   1  is   an  example   of   a   transshipment  problem  with   a  basic   feasible 
solution   specified  and  Figure   2  has    the  basis   and  the   associated  arborescence.      Our 
)ictorial  representation  with   the   "root"   at   the   top   rather   than   at   the  bottom  is 
fairly  standard. 

For  node      i,      p(i)      is    called   the  $ft.<l&<LQJLb>t>OK.  of     i      and   the   rooted  arborescence 
jLs    called  the  phJLdfLCJLkkOh.  QHa.ph.      The   predecessor  graph   is    closely   related   to  but   not 
Identical  with    the   classical  result    of  Koopmans    [40]    that    the   arcs   of   a   transshipment 
)asis    form  a   tree   over   the  nodes   of   the   problem.      The    classical   tree  preserves    the 
orientation   of   arcs   in    the   original  network  and  does  not    include  node     m  +   1. 

The  predecessor   function      p(    )      is   a  well  known   compact  way   to   represent   trees 
ind  rooted  arborescences    and  has  been   used  in  network   algorithms    for   at   least    15  years, 
.'he  predecessor  graph  has   often  been   used  interchangeably  with   the   classical  basis 

ree.      It   is   important    to   distinguish  between   them  because   the   predecessor  graph   is   a 
lata  structure   that  supports    the   computation   of   the   algorithm  and  the   orientation   of 

he   arcs   indicates   the   unique   direction   to    the   root   rather  than   a  direction   in   the 
tetwork.      Furthermore,    the   predecessor  graph   can  be   extended   to   triangular  bases    that 
Lave  no   underlying  network. 
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Figure  1.   A  Single  Commodity  Transshipment  Problem. 
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Figure    2.      A   Transshipment   Basis    (for   Figure    1) 
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The   predecessor  function     p(   )      can  be   iterated  to   construct   the  unique  path 
from  any  node    to   the   root;   this   path  is   called  a  backpaXk.      The  immzdicvtQ,  6u.C.deAi>0Kt> 
of   a  node,    if   any,    are   the   first  nodes   encountered  on   all  paths  except   the  backpath 
to   the   root   and  all  the  nodes   on  these  paths   are   called  the  i> UL£(lQJ>i> OHA .      Another 
characterization   of  the   successors   of  node      i      is   that   they   are   all  nodes  whose 
backpath   to   the   root   includes  node      i.      A  tree   can   also  be   represented  by  the 
immediate   successors   of  each  node;   however,    since   the   class   of  trees    that   arise   in 
network  problems   is    called     m-ary     with     0      to     m     immediate   successors    for  each  node, 
this   is  more   difficult    to  maintain   dynamically   than   the   predecessor   function  which 
always  has   a  single   unique  value   for  each  node  except   the   root. 

In  general,    there   are  many   different   triangulations    for  any  given   transshipment 
basis.      (Note   that   at   each  step   in   the   construction   there  may  be  several   choices   for 
the  next    column.)      However,    all  such   triangulations  yield   the   same  predecessor  func- 
tion  and  graph    (where   the   ordering  of  successors   right    to    left    for  any  node   is 
immaterial).      Thus,    the  predecessor  graph   does   not    completely   represent   a  triangula- 
tion  without    additional  information,   namely     an    ordering  of   the   rows.      A  mathematical 
development   of    these   properties    and  their  implications   is   given   in   the   companion 
paper. 

The   relationship  between   the   algebraic  view  of  the   simplex  algorithm  and  the 
graph   theoretic  view  of  much   of  the  network   literature   can  be   shown  by   describing  the 
operations   of   the   simplex  algorithm  and  the   triangulation   in   terms   of   the  predecessor 
graph.      A  graph    theoretic  proof    of   the   triangulation  theorem  identifies   a  node  with  ea 
row  and  an   arc  with  each   column  with   two  nonzero  elements.      Also   included  is   the   root 
node;    columns  with   a  single   -1    (or   1)    are   represented  as    an   arc   to    (from)    the   root 
node.      A  triangulation   and   the   predecessor  graph   are   constructed  by   first   selecting 
the   root  node.      Select   any   arc  with  one  end   the   root   and  add   the  new  node   and  arc 
orienting   the   arc   toward   the   root.      For  the     k —     step   of  the   construction     k     nodes 
and     k  -   1     arcs   are   already   in    the   graph.      Select   an   arc  with   exactly  one  node 
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already  in  the  graph,  add  the  new  node  and  the  arc  to  the  graph  orienting  the  arc 
toward  the  old  node.   The  resulting  triangulation  is  characterized  by  the  order  in 
tfhich  the  nodes  and  arcs  are  introduced. 

For  the  network  linear  program  each  node  is  identified  with  an  equality 
constraint  and  each  arc  k   from  node  i   to  node  j   is  identified  with  a  variable 
tfith  column  N   having  a  1  in  row  i,   a  -1  in  row  j   and  O's  elsewhere.   For  the 
discussion  below  it  is  assumed  that  the  basis   B   is  in  triangular  form  and  to 
simplify  notation  it  is  momentarily  assumed  that  row  i   of  B   corresponds  to  node  i 
Eor  all  i  e  W.   (Renumber  the  nodes  if  necessary  so  that   I  =  {1,2 , . . . ,m}.)   It  is 
further  assumed  that  all  the  diagonal  elements  of  B   are  1;  if  B   initially  has  a 
-1  on  the  diagonal,  the  reflection  a,     -   5t   transforms  it  to  a  1.   This  may  be 
/iewed  as  transforming  the  basis  variables  to  make  the  orientation  of  arcs  in  the 
jredecessor  graph  the  same  as  the  orientation  in  the  Koopmans  basis  tree. 

Given  the  vector  of  simplex  multipliers   u,   the  priceout  formula  (11)  for  non- 
jasic  arcs  (step  SI)  simplifies  to   c,  -  u.  +  u..   Thus  the  priceout  involves  only 
addition  operations. 

For  the  determination  of  the  arc  to  leave  the  basis  in  step  S2  the  system  of 

k    k  k 

equations   BZ  =  N   must  be  solved  for  Z  .   This  calculation  is  described  by 

showing  how  to  solve   BQ  =  e.   for  Q  .   (Q    is  thus  the   j —  column  of  the 

Inverse  of  B.)   Since  B   is  triangular,   Q   can  be  obtained  by  simple  backward 

solution:   the  m,m-l, . . . , j+1  elements  of  Q   are  seen  to  be  0,  the  j —  element 

Ls  1.   Setting  the   j —  element  equal  to  0  puts  a  1  in  the  modified  righthand  side 

Ln  the  row  corresponding  to  the  off diagonal  -1  (if  any)  in  the  j —  column  of  B, 

this  row  is   p(j).   As  before,  elements   j+1 , .  .  .  ,p(j  )-l  of  Q    are  0  and  the   p(j) — 

element  is  1.   This  continues  putting  1  in  the  p  (p  ( j) )  ,p(p  (p  ( j) ) )  ,  etc.  elements  of 

}   until  a  column  with  no  nonzero  offdiagonal  is  encountered.   Elements  back  to 

elements   ...,2,1  are  then  set  to  0. 
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Thus,    in  terms   of   the  predecessor  graph,    all  arcs   traversed  on   the  backpath  from  | 
node      j      to   the   root  have   an  element   1  in     QJ ;    all  other  elements   are  0.      Therefore 
can  be   generated  directly   from  the  precedence   function     p(   ),  which    (with     I)    gives 

the  substitution   rules   for  the  back  solution. 

k  k  i  j 

The    calculation  of     Z        follows   immediately  since     Z     =  Q     -  Q    .      The  element   in 

k  i  j 

a  row  of     Z        is    1   for  all   rows  with   a  1  in     Q        alone,    -1   for  all  rows  with   1  in     Q 

alone   and  0   for  all  rows  with  0   in  both   or  1  in  both.      Since  the  nonzero  elements   of 
Z        are   1  or  -1,    the   calculations   in  steps   S2b   and  S2c  also  involve   only  addition   and 
subtraction.      Further,    the   calculation   is    usually   reduced  enormously  by   the  extinguish- 
ment  of   the  elements    common   to  both      i      and     j      backpaths. 

As   an  example,    consider   the   triangulated  basis   in   Figure   2    (noting  that   the  nodes 

k  9  9 

have  not  been   renumbered)   with   column     N        associated  with   arc   (1,9),      BQ     =  e    ,      and 

1  1 

BQ     =  e    . 

I   =   {8      1      7      5      4      3      2   11     9      6  12   10} 

p(    )    =   {8     3     4     8     8     2     1   13     2      3     2     6} 

Q1   -    [1     1      0     0     0     0      0     0      0     0      0     0]' 

Q9=    [10001110100     0]' 

Zk  =    [0      1      0     0-1-1-1     0-1     0      0     0]' 

In  step  S3,  the  primal  solution  is  undated  using  Z  .   Also,  as  shown  by  (15), 
the  simplex  multipliers,   u,   can  be  updated  rather  than  calculated  at  each  step. 
The  algebraic  characterization  of  the  update  applies  the  I —  row  of  the  inverse  of 
B,   Q  ,   where  the  outgoing  arc  is  the  I —  column  of  B.   The  characterization  of 
Q   above  showed  that  an  element  of  Q   is  nonzero  and  equal  to  1  only  if  the 
corresponding  basic  arc  is  traversed  on  the  backpath  from  node  j   to  the  root.   It 
follows  that  Q   is  all  0's  except  for  1  in  the  I —  element  and  a  1  for  any  element 
j   such  that  the  outgoing  arc  is  on  the  backpath  from  j   to  the  root.   In  terms  of 
the  predecessor  graph,   Q.   is  1  for  node   £   and  each  of  its  successors  and  0  else- 


18 


'here.   Since  the  nonzero  elements  of   Z   are  1  and  -1,   A   in  (15)  is  either  plus 
»r  minus   cl  -  u.  +  u.   where  k  is  the  incoming  arc.   Since  Q   is  all  O's  and  l's, 
:he  update  of   u  is  accomplished  by  adding  A   to  u   and  to  the  u's  of  (only)  the 
uccessors  of  node  I. 

As  will  be  shown  in  the  next  section,  the  simple  additive  updates  of  the  primal 
olution  and  dual  multipliers  in  step  S3  are  actually  accomplished  in  a  single  inte- 
;rated  process.   This  ("pivot")  process  simultaneously  performs  the  updates  while 
■etriangulating  the  new  basis.   Fortunately,  the  retriangulation  is  also  inexpensive 
o  perform  on  a  transshipment  basis  with  a  single  new  column  violating  triangularity. 

IMPLEMENTATION 

The   design   of   large   scale  programming  codes  necessarily   involves  many  significant 
lecisions  which  have   major  impact.      The   following   fundamental  principles  were   used  to 
■esolve   design   questions   in   the   development   of   the   code   reported  here. 

The   code   is   designed   for   large   scale  problems.      Even   though  experimental 
testing  will  be   confined  by   economic   considerations    to  problems  with  some 
arbitrary  maximum  size    (say,    10,000  equations),    the   design   decisions   should 
lead  to   a   code  with  superior   large   scale  performance. 

1.      The   goal  is   a   code    for   the  most   general   capacitated  transshipment   problem. 
While   problems  with   specializations    (e.g.,    uncapacitated,    transportation, 
assignment)    must  be   solved,    the  basic   code  will  not  be   tailored  to  these 
special   features.      In   addition,   no   special  numbering  of  nodes,   extensive 
preprocessing,    or  other  design   specificity  will  be   required   that  will  limit 
the    capability   of   the    code.      Efficient   solution   of  problems   should  not 
require   detailed  advance   knowledge   of  problem  structure    (for  instance,    a 
feasible   initial  solution) .      Problems  with  multiple   arcs  will  be   accommo- 
dated. 
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3.  For  practical  problems  the  number  of  arcs  is,  in  general,  much  greater  than  the 
number  of  nodes,   m.   However,  the  problems  are  usually  sparse  in  the  sense 
that  the  number  of  arcs  seldom  approaches  the  maximum  number  m(m-l)  of  oriented 
node  pairs.   Thus,  it  is  significantly  more  expensive  to  store  information 
associated  with  each  arc  than  it  is  for  each  node  and  prohibitively  expensive 
to  store  a  node-arc  incidence  matrix.   Practical  general  minimum  cost  network 
flow  problems  are  always  heavily  capacitated. 

4.  It  is  important  to  produce  a  code  that  is  machine  independent  as  well  as 
efficient.   For  example,  machine  specific  features  such  as  assembly 
language,  use  of  particular  offline  mass  storage  devices,  storage  of  data 
using  bit  string  logical  vectors,  use  of  other  architectural  curiosities 

on  particular  machines  or  nonstandard  language  features  are  all  to  be  avoided. 
The  language  used  is  basic  FORTRAN. 

5.  Where  feasible,  speed  of  execution  will  be  given  preference  over  economy  of 
space  for  data,  storage. 

6.  Since  the  program  will  be  used  for  comparisons  of  various  data  structures  on 
a  wide  variety  of  network  problems,  it  must  be  equipped  with  effective 
external  tuning  parameters.   While  some  tuning  is  possible  for  the  pivot 
mechanism,  the  pricing  scheme  invites  especially  close  scrutiny  for  tuning 
purposes . 

Once  the  design  of  an  efficient  network  code  is  chosen,  consideration  will  also 
be  given  to  additions  of  other  advanced  features  such  as  "in-core/out-of-core" 
operation,  implicit  arc  generators,  crashed  bases,  nonlinear  costs,  post  optimal 
analysis,  and  so  forth.   These  extensions  will  not  be  allowed  to  interfere  with  the 
basic  design  goals,  but  they  should  not  be  precluded  by  the  basic  design  decisions. 
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The   description   of   GNET  begins  with  discussion   of   arc  and  node   storage 
representation.      Next,    the   ratio   test   is    described.      The  pivot   is   explained  function- 
ally   (with  a  more  extensive   algebraic  derivation   left    for  the   companion  paper) . 
Finally,    the   pricing  mechanism  is   examined.      This   order  is    chosen    (steps   S2   and  S3 
followed  by   SI)    to   faithfully   report   the   implementation  historically   and  to   lead 
smoothly   to   the   computational  performance   tests.      Hereafter,   notation  with   upper-case 
Roman   letters   indicates   a  program  variable   and  addition   of  parentheses   denotes   an 
array.      For   instance,    the   predecessor  <Wicuj  is    referred  to   as     P(   ). 

Since   there  will  be  many   arcs,    it   is    critical   to  minimize   the   stored  data 
describing  each   arc.      A  typical  input   format    (for  example,    SHARE)    for  each   arc  is: 
tail,   head,    cost,    lower  bound  and   upper  bound  on   flow.      The   lower  bounds   are   removed 
by  transformation.      If   the   arcs   are   sorted  so   that   all  arcs  with   the  same  head  are 
stored  in   contiguous   space,    the   list   of  heads    can  be   replaced  by   a  node-length   array 
whose      j —     element    is    the    location   of  the   first   arc  with  head     j.      Since  network 
models   have  many  more   arcs   than  nodes,    this   reduces   the   storage   requirements   for    the 
algorithm.      Thus,    the  network  is   stored   as    three   arc-length   arrays:      the   tails     T(    ), 
the   costs      C(   )      and  the   upper  bounds    (capacities)      CP(   );    also,    one  node-length   array 
is   used,    the   head   entries   H(    )    into   T(    ).  Positive   capacities   are   required  after 

transformation   of    lower  bounds    for   all  arcs — uncapacitated  arcs  have   capacity  set    to 
some  value   greater  than   the   total  supply.      Arcs   out   of  the  basis    at   their  upper 
bound  are  marked  with    a  sign  bit   on   the    capacity      (-CP(    )). 

It   is  natural   to   associate  with   each  node      i      (except   the   root)    the   unique  basic 
arc   that    connects      i      to  its   predecessor.      It    is    convenient   to  have   the  basic  arcs 
oriented  the   same   as   in   the  predecessor  graph,    that  is,    from  a  node   to   its   predecessor. 
This   is    accomplished  by   reflecting  arcs   as  necessary.      The  predecessor  array     P(   )      is 
marked  with   a  minus   sign   for  all  arcs   that  have   not  been   reflected.      (Subsequently 
when   using     P(    )      we  will  assume   for  simplicity  that   it   is   positive.)      The   flow  on   arc 
(i,P(i))      is   stored  in     X(i) . 
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In  step  S2  of  the  simplex  algorithm,  it  is  necessary  to  compute  capacity  minus 
flow  for  basic  arcs  with   z.  .  <  0.   It  is  convenient  to  use  a  node-length  array  to 
speed  up  this  calculation.   One  approach  is  to  store  for  each  basic  arc  a  pointer  to 
its  capacity  in  the   CP(  )   array.   Another  technique  is  to  store  the  capacity  rather 
than  the  pointer.   A  third  method  is  used  in  GNET.   The  capacity  minus  flow  is  stored 
in  a  node-length  array   CPX(  ) . 

The  simplex  multipliers  are  stored  in  a  node-length  array  U(  ) .   Figure  3  shows 
these  arrays  for  the  basis  in  Figure  2. 

After  the  incoming  arc  has  been  chosen  in  step  SI,  the  outgoing  arc  is  determine< 
in  step  S2  and  the  data  structures  are  updated  in  step  S3.   Let   k  be  the  incoming 

arc  going  from  node   i   to  node   j.   The  possible  outgoing  arcs  correspond  to  the 

k  k 

nonzero  entries  in  Z  .   We  have  already  seen  that  Z   is  zero  for  all  nodes  common 

to  both   i   and  j   backpaths.   Let  the  join   be  the  first  node  on  the  backpath  from 

i  to  the  root  that  is  on  the  backpath  from  j   to  the  root  (the  join  is  the  extin- 

guishment  point  for   Z  ) .   The  possible  outgoing  arcs  are  the  arcs  on  the  backpath 

from  i   to  the  join  and  the  arcs  on  the  backpath  from  j   to  the  join. 

It  is  critical  to  identify  the  backpaths  from  i   and  j   to  the  join  efficiently 
The  trial  and  error  of  the  classical  "stepping  stone"  methods  of  most  textbooks  will 
clearly  not  suffice  for  any  but  trivial  problems.   We  discuss  four  methods  to  identify 
the  backpaths  from  i   and  j   to  the  join  node.   Note  that  only  predecessor  infor- 
mation is  available  to  the  program  at  each  node  and  that  our  data  structure  has  no 
global  view  of  the  arborescence  as  does  the  reader  of  Figure  2. 

The  most  naive  method  is  to  mark  in  some  way  all  the  nodes  on  the  backpath  from 
i   to  the  root.   Then,  the  first  marked  node  encountered  on  the  backpath  from  j   to 
the  root  is  the  join.   This  method  is  satisfactory  for  smaller  problems,  but  for 
larger  problems  it  is  more  efficient  to  avoid  the  iteration  along  one  complete  back- 
path  all  the  way  to  the  root  by  keeping  additional  information  about  the  tree.   Also, 
in  this  way  it  will  not  be  necessary  to  unmark  the  marked  nodes  before  performing  the 
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Figure    3.      GNET  Arrays    (for   Basis    in   Figure    2). 


f 


TTwo  candidate  queue  arrays  also  used  by  GNET  are  omitted  for  clarity.   Actual 
simplex  multipliers  in  GNET  would  all  be  exactly  1,188  units  smaller  than 
shown.   This  difference  is  the  high  cost  associated  with  artificial  demand 
arc  13->8,  initially  computed  as   BM=  Nodes  x  maximum  |  cost  |  . 
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next   simplex  pivot. 

One  efficient  method  is   to  store   for  each,  node   the  number  of  nodes   on   the 
backpath   to  the   root,    call   the   d&ptk  of  the  node,      D(   ) ,      in   the   tree.      With   depth 
information   available   it   is   possible   to  iterate   the  backpaths   synchronously   from     i 
and     j      to  identify   the  join  without   iterating  either  backpath  past   the  join.      Depth 
is   used   to  indicate  which  backpath  node   is    deeper  in   the   tree   and  should  be   iterated. 
When  both  backpath  nodes  have  matching   depths,    the  nodes   are   compared   for  equality. 
A  match   indicates   the   join,    and  a  mismatch   indicates   that  both  backpaths   should  be 
iterated  for  another   comparison. 

Another  efficient  method  similar  to   the   depth   approach   is    to  store   for  each  node 
the  number   of  nodes    in   its   subtree,    called   the  nwfodh.  o{  i>u.C(L2Ai>0hA  ,      NS(    ).      Startin 
with  nodes      i      and     j,      the  backpath  node  with   strictly   fewer  successors   is   iterated. 
When  both  backpath  nodes   have   the   same  number  of  successors,    a  match   indicates   the 
join   and  a  mismatch   forces   iteration   of  both  backpaths. 

The   fourth  method  will  be   discussed  below. 

The   latter   three   methods    for   locating   the  join   look   only   at   arcs   that   are   on   the 
backpaths    to   the  join,    thus   it   is   possible   to   determine   the   outgoing  arc  while 
searching   for   the   join.      As  noted  above,    all  the   arcs    on   the  backpath    from     i      to  the 
join   are   the  +1  elements   of      Z        and  all   the   arcs    on   the  backpath   from     j      to  the 
join   are   the   -1  elements.      The   ratio   test   step   S2   is   then  simply 

CP(k)        the   capacity   of   the   incoming  arc 
min   •(  X(    )  for  arcs    on   the  backpath   from     i      to   the   join 

CPX(    )      for  arcs   on   the  backpath   from     j      to  the   join. 

The   computational  simplicity   of   this    ratio   test   is    the   rationale   for   the   reflection 
of  basic   arcs   and   the   adoption  of   CPX(    ) . 

If    the   incoming   arc  is   out   of   the  basis   at   its    capacity,    then  step   S2  may  be 
viewed  as   increasing  flow   in   a  fictitious   arc   from  node     j      to  node      i     with   the  same 
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capacity. 

A  major  proportion  of  the  work  of  each  simplex  step  is  S3,  the  pivot.   In  this 
step,  the  entering  and  leaving  arcs  are  exchanged,  the  flows,   X(  )   and  CPX(  ), 
are  updated,  the  simplex  multipliers,   U(  ),   are  changed,  and  the  simplex  data  arrays 
P(  )   and  D(  )   (using  the  depth  mechanism  for  example)  are  modified. 

Within  this  step  the  coordination  and  sequencing  of  operations  are  critical  to 
the  efficiency  of  the  network  algorithm  because  the  manipulation  of  many  nodes  and 
heavy  use  of  the  simplex  data  structure  are  involved.   If  properly  done,  this  step  is 
the  elegant  central  part  of  the  code  that  can  be  executed  by  a  computer  quickly; 
however,  the  explanation  will  be  somewhat  intricate — this  part  of  the  algorithm  is 
considered  by  many  to  be  the  "secret"  part  of  the  proprietary  codes. 

To  illustrate  a  typical  pivot,  Figure  2  shows  the  entering  arc   (i,j),   join  and 
leaving  arc   (c,d).   Call  the  backpath  from   j   to   c   the  p-ivot  itum.      Figure  4  in- 
cludes the  predecessor  graph  after  the  pivot.   Notice  that  the  subarborescence  with 
root   c  is  "rehung"  from  node   i,   and  that  nodes  on  the  pivot  stem  have  their 
predecessor  relationships  reversed.   The  flows,   X(  )   and  CPX(  ),   change  only  on 
the  backpaths  from  i   and  j   to  the  join.   The  simplex  multipliers,   U(  ),   and  the 
depth,   D(  ),   are  recomputed  for  all  nodes  in  the  subarborescence  with  root   c. 

The  most  expensive  part  of  the  pivot  is  the  update  of  the  simplex  multipliers  by 
the  addition  of  X     to  the  U(  )'s   of  node   c   and  all  nodes  in  the  subarborescence 
with  root   c.   With  the  data  structure  presented  so  far,  it  is  not  easy  to  identify 
all  the  (successor)  nodes  in  a  subarborescence.   The  identification  of  these  nodes 
can  be  facilitated  by  a  £Aa.\)<LKkaJl   data  structure  that  begins  at  the  root  of  the 
predecessor  graph  and  exhaustively  "walks"  through  all  the  nodes  in  the  same  sequence 
that  the  nodes  occur  in  a  triangulation  of  the  basis.   This  is  done  with  a  node-length 
array   IT(  )   whose   i —  element  is  the  next  node  to  visit  from  node   i.   IT(  )   is 
thus  a  different  way  to  represent  the  information  in   I.   It  is  convenient  to  make 
this  a  circular  list  by  setting  the   IT(  )   of  the  last  node  in  the  triangulation 
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Figure    4.      GNET  Arrays   and   Basis   after   Pivot    (Figure    2) 
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equal  to  m  +  1. 

From  the  triangulation  construction  it  is  clear  that  a  node  must  be  selected 
before  any  of  its  successors.   The  construction  of  the  triangulation  can  always  be 
modified  to  select  all  the  successors  (if  any)  of  a  node  before  any  other  node  is 
considered.   For  this  restricted  class  of  triangulations  the  corresponding  traversal 
IT(  )   will  visit  all  the  nodes  of  a  subarborescence  contiguously,  precisely  as  is 
required  in  the  update  of  the  simplex  multipliers.   This  traversal  is  the  obvious 
extension  to  the  m-ary  case  of  the  phJLOKdSLh.   (or,  equivalently ,  dynastic  order) 
traversal  in  computer  science  literature.   The  recursive  definition  of  the  preorder 
traversal  given  in  the  computer  science  literature  [41]  reveals  precisely  its  value 

in  updating  the  simplex  multipliers. 

Let's  look  at  the  work  that  must  be  done  in  the  pivot.   The  algorithm  visits 
each  node  of  the  subarborescence  with  root   c  exactly  once.   It  proceeds  up  the 
pivot  stem  one  node  at  a  time.   At  each  stem  node,  the  successors  of  the  next  lower 
stem  node  have  already  been  visited.   The  unvisited  successors  of  the  current  stem 
node  can  be  divided  into  two  groups:   the  nodes  visited  in  preorder  (by  iterating 
IT(  ))   from  the  stem  node  until  the  next  lower  stem  node  is  encountered,  called  the 
left  successors  of  the  stem  node,  and  the  remaining  unvisited  nodes  in  preorder  called 
the  right  successors  of  the  stem  node.   For  example,  in  Figure  2  the  left  successor  of 
stem  node  2  is  11,  and  the  right  successors  are  6  and  12.   Figure  5  gives  a  general 
description  of  the  pivot. 

We  have  experimented  with  three  different  data  arrays  that  will  support  the 
0ne-po64  pivot   computation:   depth   D(  ),  number  of  successors  NS (  ),   and  an 
additional  structure   PD(  )   to  be  discussed  below.   We  did  not  compare  another,  less 
efficient  two-pass  method  which  marks  nodes  with  sign  bits,  and  later  unmarks  them, 
much  like  the  procedure  described  for  locating  the  join  node.   These  arrays  are  used 
in  the  pivot  solely  to  answer  the  question  (Figure  5)  "An  unvisited  right  successor 
remains?";  this  question  is  nontrivial  because  the  node  information  now  available  to 
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us  is  local  in  nature. 

Figure  6  shows  a  detailed  description  of  the  pivot  traversal  using  depth.   The 
purpose  of  visiting  each  node  is  to  update  U(  ).   The  pivot  also  updates  the  arrays 
D(  )   and  IT(  ),   but  these  operations  are  omitted  for  clarity.   The  arrays  P(  ), 
X(  )   and  CPX(  )   are  recomputed  for  the  pivot  stem  only.   The  update  of  IT(  )   is 
easy  because  it  changes  only  for  the  nodes  on  the  pivot  stem  as  well  as  for  the  last 
left  successor  (if  any)  and  the  last  right  successor  (if  any)  of  each  pivot  stem  node. 
For  nodes  on  the  pivot  stem,   D(  )   is  updated  as  the  pivot  moves  up  the  stem.   The 
right  and  left  successors  of  each  pivot  stem  node  "inherit"  their  depth  from  the 
pivot  stem  node — if  the  depth  of  the  pivot  stem  node  changes  by  DADJ,   then  so  does 
the  depth  of  its  right  and  left  successors. 

The  number  of  successors,   NS  (  ),   can  also  be  used  as  illustrated  in  Figure  7. 
The  updating  of  the  number  of  successors  is  not  shown  in  the  figure.   The  number  of 
successors  is  easy  to  update  because  it  changes  only  for  the  nodes  on  the  backpaths 
from  i   and  j   to  the  join.   Thus,  the  update  can  be  performed  for  stem  nodes 
simultaneously  within  the  pivot,  and  for  the  other  backpath  nodes  (from  i   and  d) 
with  the  X(  )   and  CPX(  )   flow  update. 

Degeneracy  is  a  critical  issue  in  transshipment  problems.   In  some  of  our  test 
problems  more  than  90%  of  the  (tens  of  thousands  of)  pivots  are  degenerate.   The 
search  for  the  join  may  be  aborted  when  degeneracy  is  encountered  since  the  only 
purpose  of  the  search  is  identification  of  the  leaving  arc  and  backpaths  for  flow 
change  (zero  in  the  case  of  degeneracy) .   Stopping  short  also  tends  to  make  the  number 
of  nodes  visited  by  the  pivot  smaller.   In  the  depth  version,  the  lowest  degenerate 
arc  is  chosen  to  leave,  while  in  the  version  with  number  of  successors  the  smallest 
number  of  nodes  is  chosen  for  the  pivot.   But  the  successors  version  must  still 
locate  the  join  to  update  the  number  of  successors  on  each  backpath,  a  relatively 
easy  process,  while  the  depth  version  requires  no  further  backpath  search. 
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Figure  7.   Pivot  (Retriangulation)  Segment  Using  Number  of  Successors, 
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Another  data  array  that  can  be  used  in  place  of  depth  or  number  of  successors 
is  position  in  the  triangulation,  or  equivalently  position  in  IT(  ).  For  each  node  I 
i  define  pJi2.on.d2A.  dlbtayic.2.  PD(i)  to  be  the  row  number  in  the  triangulation  of 
equation  i  or  equivalently  to  be  the  number  of  iterations  of  IT(  )  beginning  witl 
the  root  to  get  node   i.   Set  PD(root)  =  0. 

The  search  for  the  join  is  particularly  easy  with  preorder  distance.   The 
following  proposition  gives  a  simple  construction  that  determines  the  join. 

# 

Proposition  (McBride  [44],  Graves  and  McBride  [32]):   Given  a  basis  with  a  preorder 
triangulation,  for  any  two  nodes   i   and  j   such  that  PD(j)  >  PD(i)   the  first  node 
h   on  the  backpath  from  j   with  PD(h)  ^  PD(i)   is  the  join. 

Proof :   The  construction  always  determines  a  node  since  PD(root)  =  0.   The  join  is 
the  first  node  on  the  backpath  from  one  node,  say   j,   that  is  also  on  the  backpath 
from  the  other  node   i.   Thus  we  need  to  show  that   i   is  either  h   or  a  sucessor 
of  h   and  that   i   is  not  a  successor  of  any  other  node  on  the  backpath  from  j   to 
h.   In  any  triangulation  PD(k)  >  PD(P(k))   for  any  node  k,  thus  by  construction  h 
is  the  only  node  on  the  backpath  from  j   to  h   that  could  be  the  join.   In  a 
preorder  triangulation,  the  successors  of  any  node  k  have  contiguous   PD(  )   numbei 
beginning  with  PD(k)  +  1.   Since  PD(h)  ^  PD(i)  <  PD(j) ,   i  either  equals  h  or  is 
a  successor  of  h.   Therefore,   h   is  the  join  of  i   and  j.    □ 

Using  PD(  )   the  outgoing  arc  can  be  determined  simultaneously  with  the  search 
for  the  join  and  if  degeneracy  is  discovered,  the  search  for  the  join  is  stopped. 

The  question  in  the  pivot  "An  unvisited  right  successor  remains?"   in  Figure  5 
is  equivalent  to  "PD(P(NT))  ^  PD(NR)?"   and  the  flowchart  for  preorder  distance  is 
equivalent  to  Figure  6  with   "D(NT)  >  D(NR)"   replaced.   PD(  )   is  recomputed  for  alJ 
nodes  in  the  subarborescence  with  root   c  and  some  additional  nodes. 

PRICING  MECHANISM 
The  pricing  operation  of  step  SI  in  the  simplex  algorithm  requires  a  great  deal 
of  computational  effort.   As  in  other  large  scale  mathematical  programming  problems, 
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letwork  codes  can  spend  more  than  half  of  their  execution  time  selecting  incoming 
variables  by  pricing.   Thus,  the  pricing  mechanism  is  crucial  to  overall  performance 
and  must  also  provide  the  flexible,  broad  and  effective  external  controls  necessary 
to  permit  tuning  of  the  code  for  network  problems  with  special  or  even  bizarre  struc- 
ture. 

Choice  of  pricing  strategy  is  truiy  an  art  for  large  scale  mathematical 
programming  codes.   It  is  based  on  intuition,  experience  and  empirical  experiments 
iWith  the  class  of  problems  to  be  solved.   Although  simplex  pricing  has  received  very 
little  exclusive  attention  in  the  literature,  fast  primal  codes  all  employ  some  form 
of  multiple  or  partial  pricing  of  subsets  of  the  variables  at  each  pivot  [47]. 

Examination  of  many  problem  trajectories  suggests  that  pricing  be  performed  in 
three  major  phases — an  opening  gambit,  a  middle  game  and  an  end  game.   Initially,  a 
feasible  solution  is  being  constructed  and  pricing  must  select  incoming  variables 
carefully  among  many  available  choices  with  the  view  of  satisfying  violated  constraints. 
Later,  normal  gains  toward  optimality  seldom  justify  extensive  competition  among 
candidates.   Finally,  a  favorable  candidate  becomes  quite  rare  and  thus  has  consider- 
able value  since  eventually  all  variables  must  be  exhaustively  priced  to  verify 
optimality. 

GNET  performs  pricing  by  selectively  using  a  QZYiZhaL  i>can   mechanism  and  keeps  a 
set  of  good  pivot  information  in  a  candidate.  quLCim.      For  a  given  head  node,  the  scan 
selects  the  single  incident  arc  (if  any)  pricing  most  favorably  (for  each  tail  examined 
incident  to  a  particular  head  node  only  one  addition  and  one  comparison  are  required) 
and  places  it  on  the  candidate  queue.   The  candidate  queue  is  a  varying  length  cyclic 
list  of  "interesting  node"  and  "good  arc"  information.   Each  entry  includes  a  head 
node  and  either  a  tail  pointer  specifying  the  location  of  a  good  candidate  arc  or  a 
null  pointer  indicating  an  interesting  node  which  has  not  yet  been  priced.   The 
candidate  queue  mechanism  provides  for  user  control  of  network  solutions.   Although 
there  are  many  special  uses  for  these  controls,  for  example,  basis  crashing,  it  is 
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even  more  important  to  provide  robust  automatic  procedures  for  initialization, 
selection  of  incoming  arcs,  and  queue  maintenance  for  the  most  general  class  of 
transshipment  problems. 

Since  GNET  is  designed  for  general  capacitated  transshipment  problems,  heuristi< 
for  advanced  feasible  starting  solutions  are  not  employed.   These  heuristics  are 
usually  designed  for  problem  specializations  and  the  cost  of  their  use  is  not  always 
clearly  justified  for  large  problems.   GNET  uses  an  initial  basis  of  artificial  arcs 
connected  to  the  root  node  with  initial  flows  equal  to  supplies  and  demands  and  a 
high  cost  attached  to  each  artificial  arc  incident  to  a  demand  node.   The  candidate 
queue  is  initially  loaded  with  all  these  demand  nodes. 

For  each  pivot,  the  incoming  arc  is  selected  from  the  candidate  queue  by 
examining  entries  in  a  block  of  specified  size  ("number  examined"  NNE) :   each 
candidate  arc  encountered  is  individually  repriced;  each  interesting  node  is  priced 
out  by  the  general  scan.   The  most  favorable  arc  found  in  the  block  (if  any)  is 
chosen  to  enter  the  basis.   Also,  all  arcs  pricing  favorably  are  returned  to  the 
candidate  queue. 

The  opening  gambit  is  designed  to  accelerate  the  achievement  of  feasibility. 
For  this  sequence  of  ("number  starting"  NNS)  pivots  the  high  cost  of  infeasibility 
is  likely  the  cause  for  arcs  to  price  favorably.   These  pivots  essentially  build 
chains  from  demand  nodes  to  supply  nodes.   Accordingly,  GNET  stimulates  such  early 
chain-building  by  returning  to  the  candidate  queue  the  head  and  the  tail  of  each 
incoming  arc.   Thus,  the  relatively  costly  general  scan  mechanism  is  directed 
subsequent  to  the  pivot  to  specifically  price  out  arcs  connected  to  these  two 
interesting  nodes  probably  representing  infeasible  problem  constraints.   This  "demanc 
driven"  scheme  works  particularly  well  on  problems  with  relatively  large  numbers  of 
demand  nodes  and  actual  structural  chains.   It  also  motivates  the  arc  list  organiza- 
tion by  head  node  rather  than  the  usual  tail  ordering. 
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During  the  opening  gambit,  if  no  favorable  arc  is  discovered  after  examining  a 
block  on  the  candidate  queue,  another  block  of  entries  is  accessed,  and  so  forth, 
until  an  incoming  arc  is  discovered  or  the  candidate  queue  is  completely  emptied.   An 
exhausted  queue  is  refreshed  by  directing  the  general  scan  to  a  specified  number 
("page"   IPG)  of  head  nodes. 

After  the  opening  gambit,  the  nodes  of  the  incoming  arc  are  no  longer  inserted 
in  the  candidate  queue.   After  each  complete  cycle  around  the  queue  it  is  refreshed 
by  a  general  scan  of  IPG  nodes.   Each  of  these  scans  begins  where  the  previous  scan 
ended  and  thus  moves  cyclicly  through  the  arc  list  pricing  arcs  in  contiguous  storage 
locations.   When  managed  by  these  procedures,  the  candidate  queue  generally  grows 
during  and  just  after  the   opening  gambit  and  then  shrinks  finally  to  a  few  good 
candidates.   The  end  game  is  played  by  concentrating  on  these  last  few  candidate  arcs. 

A  major  portion  of  the  experimentation  with  sample  problems  has  been  devoted  to 
tuning  the  three  pricing  parameters  NNE,   NNS ,   and  IPG  with  the  objective  of 
estimating  optimal  settings  as  simple  functions  of  solution  progress  and  easily  measured 
problem  characteristics  such  as 'numbers  of  nodes,  arcs,  supply  and  demand  nodes, 
degree  of  capacitation ,  cost  range,  and  so  forth.   Also,  the  sensitivity  of  performance 
to  problem  structure  has  been  investigated. 

High  resolution  internal  computation  timing  is  often  difficult  to  carry  out  on 
contemporary  multiprogramming  computer  systems.   Some  published  results  are  highly 
suspect  due  to  unfortunate  oversights  in  timing  design.   The  computer  timing  routines 
often  take  longer  to  execute  than  the  network  algorithm  coding  segments  under  study. 
Therefore,  in  order  to  isolate  times  of  pricing  and  pivots,  a  special  experiment  was 
designed  to  minimize  relative  timing  error.   Sample  problems  were  run  to  completion 
by  GNET  under  fixed  parameter  settings,  each  producing  a  single  gross  solution  time 
and  number  of  pivots.   Next,  histories  of  the  entering  arcs  for  the  solutions  were  fed 
to  a  specially  modified  code,  essentially  eliminating  the  pricing  mechanism;  these 
solution  times  represent  almost  exclusively  pivot  time.   Also,  nontiming  runs  were 
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made  with  another  code  with  extensive  internal  data  collection  for  detailed  perform- 
ance analysis. 

Pricing  schemes  were  tried  ranging  from  "first  negative"  (terrible)  to  "most 
negative"  (worse) .   Optimal  intermediate  settings  of  NNE   and  IPG   lead  to  a  nearly 
equal  distribution  of  time  between  pricing  and  pivoting.   As  the  length  of  the 
opening  gambit,   NNS,   is  increased,  solution  time  is  greatly  reduced.   However, 
beyond  some  point,  times  again  go  up  as  the  candidate  queue  becomes  clogged  with  bad 
nodes  inserted  by  the  pivots.   The  size  of  the  candidate  queue  swells,  and  then 
declines  as  predicted,  with  the  maximum  size  a  complicated  function  of  all  three 
parameters.   As  expected,  several  complete  cyclic  sweeps  of  the  entire  arc  list  in 
blocks  of   IPG  head  entries  always  occur  late  in  the  end  game,  but  the  candidate 
queue  focuses  attention  on  the  best  remaining  candidates.   The  total  number  of  pivots 
is  very  sensitive  to  NNE.   However,  total  solution  time  is  relatively  stable  as  NNE 
and  the  other  tuning  parameters  are  varied  from  reasonably  good  default  settings. 
Experiments  to  change  NNE   dynamically  during  solution  progress  have  not  improved 
performance.   In  fact,  many  dynamic  tuning  schemes  have  been  tried  without  much 
success . 

Hundreds  of  such  calibration  runs  have  been  performed  in  an  attempt  to 
determine  empirically  the  form  of  the  response  surface  of  execution  time  as  a  func- 
tion of  the  pricing  parameters.   To  date,  the  best  general  automatic  default  settings 
for  transshipment  problems  are 

NNE  =  32    entries  examined  in  each  candidate  queue  block, 

NNS  =  3m/4  pivots  in  the  opening  gambit , 

IPG  =  ra/10  head  nodes  priced  to  refresh  the  candidate  queue 
in  each  general  scan  block. 

These  suggested  settings  are  surprisingly  robust  for  a  wide  variety  of  problems.  How- 
ever, for  particular  classes  of  problems  sharing  uncommon  structure,  specific  tunin 
can  achieve  remarkable  further  improvements. 
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Degeneracy  is  common  in  many  problems  (commonly  90  percent  of  the  pivots) ;  since 
degenerate  pivots  require  less  work  in  steps  S2  and  S3,  they  can  actually  be  a  bonus 
rather  than  a  worry.   The  pricing  mechanism  helps  obviate  the  need  for  any  additional 
uethods  to  avoid  cycling  by  constantly  shuffling  the  cyclic  pricing  agenda.   (No  cycl- 
ing has  ever  been  encountered  in  our  experiments.)   Should  terminal  cycling  occur, 
the  exact  integer  solutions  will  certainly  lead  to  an  infinite  computational  loop. 
In  such  an  (unlikely)  case,  specification  of  slightly  modified  pricing  parameters 
(especially  NNE)  will  almost  certainly  avoid  the  cycle.   Formal  techniques  for  deal- 
ing efficiently  with  degeneracy  and  cycling  of  primal  simplex  network  codes  invite 
further  research. 

COMPUTATIONAL  EXPERIENCE 

The  family  of  GNET  codes  first  presented  at  the  Spring  1975  ORSA/TIMS  meeting 
represents  the  state  of  the  art  in  fast  large  scale  minimum  cost  network  flow  systems. 
Tables  1  and  6  report  solution  times  and  numbers  of  pivots  for  a  set  of  standard  test 
problems  [using  NETGEN,  39]  that  have  also  been  solved  by  other  contemporary  codes. 
These  solution  times  are  achieved  with  the  pricing  parameters  set  at  their  default 
values  as  in  the  code  that  is  distributed.   Thus,  our  experiments  are  completely 
reproducible  by  other  researchers.   Benchmarks  on  various  machines  (see  Table  2,  for 
instance)  generally  agree  with  standard  hardware  comparisons  of  computer  performance 
and  show  that  the  times  in  Tables  1  and  6  are  superior  to  all  published  times  of  which 
we  are  aware.   However,  we  cannot  verify  published  claims  of  machine,  compiler  and 
installation  performance  variations.   In  any  case,  future  codes  produced  by  incorporating 
fresh  research  ideas  will  undoubtedly  make  our  current  performance  records  obsolete. 
We  do  not  hold  that  primal  network  codes  are  anywhere  near  an  asymptotic  efficiency 
level. 

The  GNET  family  presently  includes  GNET/Depth,  GNET/Successors  and  GNET /Pre order- 
distance  as  described  above,  as  well  as  other  variants  to  be  discussed  below.   The 
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TABLE  1 

GNET/Depth  Performance  on  Several  Transshipment  Examples 

(Default  Tuning) 


Problem   Nodes 


Sources    Sinks 


Arcs 


Percent   IBM  360/67 
Cap'd     Seconds 


Pivots 


NG2  7 

400 

4 

12 

2,676 

80 

2.6 

607 

NG36 

8,000 

200 

1,000 

15,000 

0 

212 

13,012 

NG37 

5,000 

150 

800 

23,000 

0 

138 

11,610 

NG38 

3,000 

125 

500 

35,000 

0 

97 

10,637 

NG39 

5,000 

180 

700 

15,000 

0.7 

113 

9,553 

NG40 

3,000 

100 

300 

23,000 

0.7 

67 

6,409 

NPS 

10,000 

50 

5,000 

21,000 

100 

441 

22,153 

XN1 

5,000 

100 

4,800 

40,000 

100 

290 

12,111 

TABLE  2 
GNET/Depth  Performance  on  'NG2  7'  -  Calibration  with  a 
Highly  Capacitated  Transshipment  Problem 


Machine 

Solut 

ion  Seconds 

CDC  7600 

0.3 

IBM  370/168 

0.5 

CDC  6600 

0.6 

TI  ASC 

0.7 

IBM  360/91 

0.8 

UNI  1108 

2.2 

IBM  360/67 

2.3 
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three  basic  versions  of  GNET  have  been  tested  on  a  suite  of  randomly  generated  problems 
and  real  formulations  to  determine  which  version  is  best.   Experiments  on  smaller  pro- 
blems (less  than  500  nodes  and  5,000  arcs)  show  that  all  three  are  remarkably  close 
with  the  depth  version  a  narrow  favorite  for  sheer  speed.   Although  preorder  distance 
is  preferable  for  mathematical  reasons  when  extending  the  data  structures  to  non- 
network  problems,  large  scale  testing  has  been  limited  to  depth  and  number  of  succesors 
versions.   GNET/Successors  execution  times  differ  from  the  GNET/Depth  times  by  at  most 
2  percent  on  the  problems  in  Table  1. 

All  versions  begin  with  a  read  and  edit  routine  and  an  arc  list  sort  to  create 
the  compressed  arc  arrays.   Supplies  and  demands  are  determined  for  each  node  and  the 
candidate  queue  is  loaded  with  the  demand  nodes.   Simplex  pricing  via  the  candidate 
queue  mechanism  is  used  to  reduce  artificial  flows  to  zero  rather  than  a  Phase  I- 
Phase  II  approach  (used  in  earlier  GNET  versions) .   The  high  cost  for  artificial  arcs 
is  computed  as  the  product  of  the  number  of  nodes  (maximum  path  length)  and  the 
maximum  absolute  arc  cost,  thus  guaranteeing  that  a  feasible  problem  will  have  no  flow 
on  artificial  arcs  in  a  final  basis.   The  cost  is  attached  only  to  artificial  arcs 
incident  to  demand  nodes;  experiments  attaching  the  cost  to  arcs  incident  to  supply 
nodes  alone,  and  to  both  supplies  and  demands  have  not  improved  performance.   GNET  is 
tuned  for  node  numbering  ascending  from  supply  through  transshipment  to  demand  nodes, 
and  departures  may  somewhat  degrade  performance  of  the  demand  driven  pricing  mechanism 
(for  general  transshipment  problems  the  general  scan  blocks  are  contiguous  and  exhaus- 
tive).  Another  underlying  assumption  has  been  that  problems  will  normally  have  many 
more  sinks  than  sources.   Problems  with  more  sources  than  sinks  or  erratic  node 
numbering  can  usually  be  easily  reformulated  to  our  preference  if  solution  speed  is 
of  prime  importance;  with  large  scale  models  this  is  a  minor  undertaking  in  the  pro- 
blem generator  software. 

Use  of  random  test  problems  in  tuning  network  codes  has  its  pitfalls.   GNET 
solves  real  network  models  much  faster  than  random  NETGEN  problems  of  nominally 
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comparable  size  and  structure;  this  suggests  that  much  remains  to  be  learned  from 
further  investigation  of  special  problem  structure.   Large  random  test  problems  are 
very  expensive  to  generate  with  NETGEN,  requiring  about  five  times  more  computer  time 
than  the  associated  GNET  optimization  (but  much  less  region).   However,  the  reproduci- 
bility of  the  experiments  and  results  is  so  important  that  it  justifies  the  expense. 
Also,  the  cost  and  awkwardness  of  using  magnetic  tapes  with  standard  problem 
libraries  is  avoided.   Unfortunately,  large  NETGEN  problems  can  vary  slightly  between 
computers  with  the  length  of  the  real  mantissa  (this  is  due  to  a  random  number  gene-  I 
rator  which  simulates  35  bit  integer  arithmetic,  normalization  of  the  integer  result 
to  a  real  0-1  variable  on  the  host  machine  and  subsequent  transformation  back  to  a 
uniform  integer  with  the  desired  range) .   NETGEN  puts  capacities  of  +1  on  all  arcs  in 
assignment  problems,  needlessly  complicating  their  solution  by  a  general  capacitated 
network  code.   Negative  demands  can  also  be  generated.   These  and  other  minor  NETGEN 
problems  are  overcome  by  a  few  program  modifications.   It  also  would  be  worthwhile  to 
add  the  capability  to  generate  multiple  arcs  and  special  problem  structures. 

GNET  execution  times  do  not  show  much  sensitivity  to  cost  ranges  contrary  to 
past  reports  [55].   Experiments  were  performed  in  which  individual  problems  were 
solved  repeatedly  with  only  the  costs  modified  either  by  low  order  digit  truncation 
or  addition  of  uniform  random  low  order  digits  as  necessary  to  provide  the  desired 
cost  range.   Solution  times  are  very  insensitive  to  cost  ranges  thus  produced,  seldom 
varying  by  more  than  15  percent.   The  outcome  of  one  such  experiment  is  shown  in 
Table  3.   This  surprising  result  appears  to  be  attributable  to  the  candidate  queue 
pricing  mechanism. 

Un capacitated  transshipment  problems  are  input  to  GNET  with  upper  bounds  on  each 
arc  exceeding  total  problem  supplies.   Minor  modifications  of  GNET  can  reduce  solutior 
times  by  about  10  percent  for  strictly  uncapacitated  problems;  a  more  important  issue 
is  the  potential  elimination  of  the  arc-length  array  of  capacities.   This  space 
economy  may  also  be  realized  on  lightly  capacitated  problems  by  modification  of  the 
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TABLE  3 


GNET/Depth  Performance  for  'NG40'  with  Varied  Ranges  for 
Random  Uniform  Costs  (Original  Range  1-100) 


IBM 

360/67 

Cost  Range 

Seconds 

Pivots 

1-10 

84 

6,726* 

1-100 

76 

6,933 

1-1,000 

74 

6,921T 

1-10,000 

74 

6.958f 

f 


Low  order  digit  truncated  from  original  problem  costs. 

Low  order  digits  generated  by  RANDU  and  concatenated  with 
original  problem  costs. 


TABLE  4 

Lightly  Capacitated  Transshipment  Problems  - 
Performance  of  GNET/Depth  and  Modifications 
Solving  the  Equivalent,  Reformulated  Uncapacitated  Problems  [11] 


NT39 


NT40 


IBM  360/67 

IBM 

360/67 

Seconds 

Pivots 

Seconds 

Pivots 

GNET (depth) 

113 

9,553 

67 

6,409 

Preprocess 

113 

9,615 

71 

7,511 

Transform  on- 

-the- 

-fly 

114 

9,652 

65 

6,733 
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capacity  array  or  by  problem  reformulation.   One  section  of  a  report  by  Cheong  [11] 
includes  a  catalog  of  problem  transformations  useful  in  network  models  and  gives 
results  of  experiments  with  several  large,  lightly  capacitated  transshipment  problems. 
Two  special  versions  of  GNET/Depth  were  prepared  which  utilize  a  well-known  transfor- 
mation to  replace  each  capacitated  arc  by  a  pair  of  uncapacitated  arcs  and  a  new  node. 
(The  tail  of  the  replaced  arc  is  moved  to  the  new  node  and  supply  equal  to  the  replace 
capacity  is  also  shifted.   The  new  node  is  then  connected  to  the  old  tail  node  by  an 
artificial  arc  with  zero  cost.)   One  modification  performs  the  transformation  to  the 
arc  list  before  solution.   The  other  carries  out  the  transformation  "on-the-f ly"  as 
arcs  with  capacities  are  introduced  into  the  basis.   Table  4  gives  an  example  of  per- 
formance for  the  codes,  with  the  modifications  each  using  one  third  less  space  for 
the  arc  arrays  with  little  speed  degradation.   As  the  proportion  of  capacitated  arcs 
increases,  the  space-time  tradeoff  becomes  much  less  favorable. 

POSTOPTIMAL  ANALYSIS,  REOPTIMIZATION  AND  FURTHER  REFINEMENTS 

In  some  applications,  analysis  of  the  sensitivity  of  the  optimal  solution  to 
modifications  in  the  supplies,  demands  and  cost  coefficients  is  desirable.   All  the 
postoptimal  analysis  for  linear  programming  problems  can  be  done  with  primal  trans- 
shipment codes.   The  data  structures  described  above  support  efficient  techniques  for 
this  analysis.   Ranging  [47]  of  problem  coefficients  traditionally  requires  informa- 
tion from  the  inverse  of  the  optimal  basis;  columns  are  required  for  ranging  of 
supplies  and  demands,  rows  are  required  for  ranging  of  costs  of  basic  arcs.   As 
described  above,  the  predecessor  function  can  generate  columns  of  the  inverse  and 
the  traversal  mechanism  can  generate  rows.   Note  that  while  the  ranging  of  cost 
coefficients  of  basic  arcs  is  simple  in  principle  and  much  faster  than  for  linear 
programming  systems,  for  applications  with  many  arcs  it  can  be  time  consuming  rela- 
tive to  the  time  to  construct  the  optimal  solution. 
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A  more  important  issue,  especially  for  applications  that  have  a  network  embedded 
in  a  larger  problem,  is  the  reoptimization  of  a  problem  after  modifications  to  the 
problem  coefficients.   The  case  for  primal-dual  or  dual  algorithms  is  often  based  on 
the  reputed  ease  of  reoptimization;  although  reoptimization  is  conceptually  easier 
than  for  primal  algorithms,  this  does  not  imply  that  primal-dual  or  dual  algorithms 
reoptimize  more  quickly.   There  have  not  been  comparisons  of  primal  and  nonprimal 

algorithms  for  reoptimization  of  large  scale  problems  (indeed,  it  is  not  clear  what 
types  of  problems  would  constitute  a  fair  comparison).   In  any  case,  the  computer 
storage  advantage  of  primal  algorithms  is  still  maintained  since  reoptimization 
requires  no  additional  data  arrays. 

The  following  design  considerations  were  used  for  a  GNET  reoptimization  procedure 
1)  there  would  be  no  modifications  to  the  primal  optimization  segments  of  GNET,  2)  no 
new  arcs  would  be  inserted  into  the  arc  list,  and  3)  the  initial  basis  for  the  reopti- 
mization would  include  as  much  of  the  previous  optimal  basis  as  possible.   Algebraic- 
ally, the  changes  to  the  supplies,  demands,  lower  bounds  and  upper  bounds  are 
translated  into  a  vector  of  changes,  d,  to  be  added  to  the  original  right  hand  side 
vector  b.   Using  the  previous  optimal  basis   B,   the  set  of  equations   Bx  =  b  +  d  is 
solved  with  artificial  arcs  introduced  into  the  basis  as  necessary  to  avoid  an 
infeasible   x.   Since  the  previous  optimal  solution  x  is  available  (b   is  not), 
Bx  =  d   is  solved  and  then  added  to   x  to  get   x. 

The  vector   d  originates  from  coefficient  modifications  by:   1)  changes  in 

supplies  and  demands,  2)  a  change  of  A   in  the  upper  bound  of  an  arc   (i,j)   out  of 

the  basis  at  this  upper  bound  (A   is  added  to   d.   and  subtracted  from  d.),  and  3)  a 

change  of  A   in  the  lower  bound  of  an  arc   (i,j)   out  of  the  basis  at  zero  (A   is 

subtracted  from  d.   and  added  to   d.).   A  traversal  mechanism  for  back  substitution 

i  3 

of   the  basis   is   provided  by   a  single   pass    through   the      IT(    )      array.      The  vector     x 

is    constructed  by  backsolving  with   the  predecessor  array     P(    ) ,      at  each   step     x.      is 

computed  as      x.    +  x. .      If      x.      is  nonnegative   and   less   than   or  equal   to  its    capacity, 
ill 
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the  backsolving  continues.   Otherwise,  an  artificial  vector  from  the  node  to  the  root 
is  introduced  to  replace  the  arc  which  leaves  the  basis  at  zero  (if  x.  <  0)  or  its 
capacity  (if  x  >   capacity) .   The  artificial  arc  is  oriented  so  that  its  flow  is 
positive  and  it  is  assigned  the  large  cost.   The  exchange  of  an  arc  and  an  artificial 
arc  involves  "rehanging"  the  node  and  all  its  successors;  this  is  accomplished  with 
only  three  changes  in  the  preorder  traversal  IT(  )  by  inserting  these  nodes  at  the 
end  of  IT(  ).   Arcs  removed  from  the  basis  are  placed  in  the  candidate  queue.  The 
backsolving  continues  until  the  whole  basis  is  recomputed.   If  any  artificial  arcs 
have  been  introduced  it  is  necessary  to  recompute  the  simplex  multipliers.   The  pro- 
blem is  then  reoptimized  by  GNET. 

If  cost  coefficient  changes  involve  only  nonbasic  arcs,  reoptimization  is 
necessary  only  if  one  or  more  now  price  favorably  (such  arcs  are  added  to  the  candi- 
date queue) .   If  cost  coefficients  change  for  basic  arcs  it  is  necessary  to  recompute 
the  simplex  multipliers  and  then  reoptimize  with  GNET. 

Since  GNET  is  so  fast ,  experience  indicates  that  if  there  are  many  coefficient 
changes  it  may  be  more  efficient  to  begin  with  an  all  artificial  basis  with  the 
previous  optimal  basis  arcs  preloaded  in  the  candidate  queue.   On  the  other  hand, 
applications  often  require  many  particular,  recurring  re optimizations  of  some  special 
type  which  permit  more  efficient  (and  problem  specific)  methods  to  be  applied. 

One  of  the  major  ideas  of  mathematical  programming  is  that  elements  of  the 
simplex  tableau  may  be  generated  as  needed  rather  than  stored  and  updated  at  each 
pivot.   For  example,  as  described  above,  in  GNET  the  columns  of  the  basis  inverse  are 
generated  when  needed  to  determine  the  outgoing  arc  by  iterating  the  predecessor 
function.  However,  the  basic  flows  and  simplex  multipliers  are  explicitly  stored  and 
the  ones  that  change  are  updated  during  each  pivot.   Analysis  of  computational  experi- 
ments show  that  the  major  portion  of  the  calculations  in  the  pivot  (step  S3)  is  to 
update  the  simplex  multipliers  U(  ) ,  depth  D(  )   and  the  preorder  traversal  array 
IT  (  )   (note  that   IT(  )   is  maintained  solely  to  allow  efficient  update  of  the 
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multipliers) .   The  authors  were  thus  led  to  consider  modifications  to  GNET  that  would 
allow  some  or  all  of  the  multipliers  to  be  generated  as  needed. 

It  follows  immediately  from  (11)  that  for  basis  arc  k   from  node  i   to  j : 
C(k)  -  U(i)  +  U(j)  =  0.   Arbitrarily  setting  U(root)  =  0,  it  is  possible  to  solve 
for  all  the   U(  )'s  by  front  substitution  with  the  triangulated  basis;  that  is, 
solve  for  the   U(  )  in  the  order  IT(root) ,IT(IT(root) ) ,  etc.   To  compute  a  particular 
U(h)   it  is  necessary  only  to  calculate  some  of  the  other  U(  )'s,   namely  those  for 
the  nodes  on  the  backpath  from  h  to  the  root.   If  k  is  the  arc  joining  node  h  to 
its  predecessor  P(h)  ,   then  U(h)  =  U(P(h))  ±  C(k)   with   "+"   if  the  arc  is  oriented 
from  h   to  P(h)   and  "-"   otherwise.   Since  there  are  many  nodes  in  large  problems 
and  only  a  relatively  few  multipliers  change  from  pivot  to  pivot,  it  does  not  seem 
worthwhile  to  generate  all  multipliers  at  each  pivot.   Rather,  we  choose  to  store 
explicitly  enough  of  the  multipliers  so  that  for  each  node   i  either  U(i)  is  stored 
explicitly  or  U(P(i))   is  stored  explicitly.   In  the  latter  case  a  single  addition 
generates  the  multiplier  only  when  it  is  needed. 

The  impetus  for  this  approach  comes  from  consideration  of  the  capacitated  trans- 
portation problem  with  many  demand  nodes  (sinks)  and  relatively  few  supply  nodes 
(sources) .   This  is  an  important  special  case  of  the  minimum  cost  network  model  that 
has  many  applications.   A  typical  model  is  a  distribution  system  with  a  few  plants  and 
many  warehouses  or  a  few  centralized  warehouses  and  many  customers.   This  is  the  type 
of  problem  most  often  encountered  in  multicommodity  distribution  problems,  e.g.  [22]. 
Another  application  that  requires  the  repeated  solution  of  many  such  problems  is  the 
traffic  assignment  problem. 

A  special  version  of  GNET/Depth  called  TNET  was  developed  for  this  problem.   The 
multipliers  are  explicitly  kept  only  for  the  relatively  few  sources  of  the  problem. 
Since  in  the  (bipartite)  transportation  problem  arcs  only  join  sources  to  sinks,  the 
predecessor  of  each  sink  node  must  be  a  source  and  thus  the  multiplier  for  each  sink 
is  computed  by  a  single  addition  (also,  the  pricing  mechanism  loads  only  sinks  as 
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interesting  nodes  during  the  opening  gambit) .   The  GNET  arc  storage  is  ideal  for  this 
calculation  since  all  the  arcs  with  sink  node  j   are  stored  contiguously;  thus  in  th 
general  scan  all  are  priced  out  together  and  only  a  single  calculation  of  the  sink 
multiplier  is  required. 

Since  the  traversal  array   IT(  )   is  maintained  solely  for  the  purpose  of 
updating  the  multipliers,   IT(  )   and  D(  )   need  to  be  maintained  only  for  the 
relatively  few  source  nodes.   It  is  also  convenient  to  include  in   IT(  )   the  rela- 
tively few  sinks  that  have  a  successor.   An  alternate  description  is  that  IT(  )   is 
maintained  only  for  the  smallest  possible  subarborescence  that  includes  all  the 
sources;  this  is  easy  to  maintain  since  at  each  pivot  the  only  possible  nodes  that 
can  join  the  subarborescence  are  the  two  ends  of  the  incoming  arc  and  the  only 
possible  nodes  to  leave  are  the  two  ends  of  the  outgoing  arc.   No  additional  storage 
is  required,  the  U(i)   for  sink  i   is  used  to  store  the  cost  of  the  arc  from  P(i) 
to   i   and  the  sinks  not  in  the  subarborescence  are  marked  with  0  in  D(  )  .   The 
predecessor  array  P(  )   is  maintained  exactly  as  before  and  the  determination  of  the 
backpaths  and  outgoing  arc  is  the  same. 

Experiments  (see  Table  5)  show  that  TNET  is  significantly  faster  than 
GNET  for  transportation  problems  with  many  more  sinks  than  sources.    The  pivot 
choice  is  the  same  in  both  the  original  and  modified  versions  so  with  the  same 

pivot  choice  is  the  same  in  both  the  original  and  modified  versions  so  with  the  same 
value  of  the  pricing  parameters  the  same  sequence  of  pivots  is  generated.   For 
NNE  =  32,   the  reduction  in  time  for  TNET  is  due  solely  to  savings  in  the  update  of 
U(  ) ,   IT(  ) ,   and  D(  )   in  step  S3  (step  S2  is  identical  and  step  SI  is  slightly 
slower  for  TNET)  .   Experiments  with  six  smaller  problems  show  that  GNET  spends  roughly 
half  its  time  in  step  SI;  with  TNET  there  was  80%  less  time  for  S3  and  5%  more  time  ii 
SI.   Since  the  pivot  step  S3  is  so  much  faster  in  the  modification  it  is  not  worth- 
while to  select  the  incoming  arc  as  carefully  as  in  GNET;  experiments  show  that 
NNE  =  8  is  a  good  setting  for  TNET.   Table  5  shows  that  with  NNE  =  8   there  are 
significantly  more  pivots  than  with  NNE  =  32,   but  total  solution  time  is  reduced. 


TABLE  5 


Transportation  Problems  with  Relatively  Few  Sources 


IBM  360/67  Seconds 
Problem   Sources   Sinks    Arcs    NNE   TNET  GNET    Pivots 


TN7 


100     4,900    40,000    32     127 

8     99 


274     11,909 
17,583 


TN8 


250     4,750    40,000 


32 

142 

285 

12,910 

8 

135 

19,759 

TN9 


10     4,990    40,000 


32 
8 


104 
79 


261 


9,574 
12,327 


Generated  with  NETGEN  [39],  all  arcs  are  capacitated, 
cost  ranges  are  0-100  (TN7)  and  0-1000  (TN8,TN9). 

TABLE  6 
XNET/Depth  Performance  with  NNE  =  8 


IBM  360/67  IBM  360/67 

Problem    Seconds      Pivots    Problem    Seconds      Pivots 


NG27 


2.7 


964 


NPS 


265 


29,045 


NG36 

111 

17,993 

XN1 

136 

19,726 

NG37 

110 

18,195 

NG38 

94 

13,124 

TN7 

112 

17,583 

NG39 

76 

14,809 

TN8 

148 

19,759 

NG40 

59 

11,002 

TN9 

92 

12,327 
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Although  TNET  is    designed   for    transportation  problems  with  many  sinks   and 
relatively   few  sources,    it   does  well    (about   the  same  as   GNET)    on  problems  with  an 
equal  number  of  sources   and  sinks.      TNET  handles   problems  with  many  sources   and  few 
sinks  by  a  minor  modification   in   the   input   that   reverses    the   orientation   of   the   arcs. 

The  same  idea  is  extended  to  the  general  minimum  cost  network  problem  in  a  pro- 
gram called  XNET.  The  multipliers  are  explicitly  kept  only  for  nodes  that  have 
successors  in  the  predecessor  graph.  For  example,  in  Figure  2  IT(  )  is  maintained 
only  for  nodes  13,  8,  1,  4,  3,  2,  6  in  that  order.  The  U(  )  and  D(  )  arrays  for 
nodes  with  no  successors  are  used  as  in  TNET.  In  order  to  keep  the  multipliers  only 
for  the  nodes  with  successors ,  it  is  necessary  to  maintain  an  array  that  records  for 
a  node  the  number  of  its  successors  that  do  not  have  explicit  multipliers  ("aggregated 
successors"      A(    )). 

The   results   in   Table   6   indicate    that  XNET   is    faster   than  GNET  on   all  problems 
and  significantly   faster  on  problems  with   relatively  many   sinks.      XNET   is   only 
slightly   slower   than   TNET  on    the  problems    that  TNET  is    tailored   for,   but   the   loss   on 
these  problems   is  balanced  by   its   generality   and  by   its   dominance   of  GNET.      XNET   is 
successful  because    the   predecessor   graphs  have  many  nodes  with  no  successors.      In  many 
practical  problems  known    to    the   authors,   most  nodes   are  pure   sinks;   since    the 
successor   in    the  predecessor  graph  of   a  pure   sink  must  be   a  node   that  is  not  a  pure 
sink,    many  pure   sink  nodes  have  no  successors. 

TNET  and  XNET  are    refinements    of   GNET/Depth;   number  of  successors   or  preorder 
distance    could   also  be   used   for  such  modifications. 

A  modified  version   of  XNET   is    also  being  used   on   a  set   of   large,    100-percent 
dense   problems.      The  network  models    are   uncapacitated  single   commodity   transportation 
problems   embedded   in   a   recent   implementation   of   a  multi commodity,    multiple    time   period 
econometric  model  described  in    [34].      A  prototype   problem  size   is   200   sources   and   300 
sinks    (and   thus   60,000   arcs)   with   a  matrix  of   region    to   region  bulk   transport   costs. 
The  XNET  modification   is   stripped  of   arc-length   arrays,   list   references   are  modified 
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for   the   cost  matrix  and  several   capacitated   features   removed.      Also,    the   candidate 
queue  is   modified   to   access   only  sink  nodes   and   to  price   a   restricted  menu  of   the 
few   (ultimately   5)    cheapest   sources    for  each  sink  during   the   opening  gambit.      Optimi- 
zation  from  a   cold  start   requires      8.6   seconds    (2350  pivots)    on   an   IBM  370/168. 
Tuning  of  pricing  parameters   produces      surprisingly      little   improvement.      Reoptimiza- 
tion  procedures   are  employed   to  exploit  period-to-period  similarity  of  optimal  bases 
(despite   significant    temporal  variations    in  problem  structure) .      A  typical   reoptimi- 
zation   time   from  a   crashed  basis    is   1.2   seconds. 

HISTORICAL  PERSPECTIVE 

Although  the  authors  have  not  had  access  to  any  other  primal  transshipment 
computer  programs,  it  is  possible  to  identify  some  of  the  major  ideas  from  papers  and 
presentations  [10,  26,  27,  28,  29,  32,  45,  46,  54,  55].   The  major  design  decisions 
and  important  coding  specifics  may  vary  widely  in  these  systems,  but  all  the  success- 
ful contemporary  large  scale  codes  seem  to  be  based  on  a  few  key  ideas.   Notwith- 
standing our  limited  knowledge  of  specific  techniques  used  by  others,  we  trace  the 
development  of  these  ideas  as  best  we  can. 

The  major  ideas  underlying  all  contemporary  primal  network  codes  are  the 
representation  of  the  basis  as  a  predecessor  graph,  a  traversal  mechanism  to  update 
simplex  multipliers,  the  use  of  depth  or  number  of  successors  or  preorder  distance  to 
determine  the  outgoing  arc  and  to  facilitate  the  update  of  the  simplex  multipliers, 
the  use  of  a  pricing  mechanism  such  as  the  candidate  queue,  and  generation  of 
simplex  multipliers. 

The  predecessor  array  is  used  in  all  the  codes.   This  construction  in  primal 
network  codes  goes  back  at  least  to  Glicksman,  L.  Johnson  and  Eselson  [25]  1960  and 
goes  back  further  in  other  disciplines.   It  is  a  standard  approach  in  many  research 
areas . 
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Srinivasan  and  Thompson   154]    1972  have  proposed   the   use   of   depth  to  identify   the 
backpaths   and  to  determine   the  outgoing  arc.      This  has  been  adopted  by   Glover,   Karney 
and  Klingman    [26]    1974,   Mulvey   [45]    1974,    [46]    1975   and  GNET/Depth  [6]    1975.      The 
method    [54]    for  identifying   the   join   of   the  backpaths    is   exactly   as   described   for 
GNET,   however,    the   update   of   depth   at  each  pivot   in   that  paper  is   a  more   involved 
approach   that   does  not   aid  in   the   pivot   as   described  here. 

Preorder   traversal  was    first   used  in    transshipment   codes   by  Glover,   Klingman   and 
Stutz    [29]    1974  who   called  it   the  Augmented  Threaded   Index  method   (ATI).      They  show 
that   it   is  more  efficient   than    the   triple   labeling  scheme   of  E.    Johnson    [37].      Indepen- 
dently,  McBride    [44]    1973   and  Graves   and  McBride    [32]    1973  have   developed   the   zero   to 
the   right    (ZTR)    triangulation   of  network  bases  which  has   since  been  shown   to  be 
equivalent   to  preorder.      Preorder  has  been   used  in   the   computer  science    literature 
for  at   least  20  years    and  has  been   used   in   contemporary  shortest  path   algorithms    [24] 
1973.      Preorder  has  been   adopted  by   Srinivasan   and  Thompson   and  Mulvey   and  it   is   used 
in   all  versions   of  GNET.      As   discussed  above   and  shown   in   Figure   6,    depth  makes 
possible   a   one-pass   update   of   the   simplex  multipliers    that  simultaneously  updates   the 
preorder   and   depth    (or  number  of  successors   or  preorder  distance) .      Since   the  pivot 
takes   much  more   time   than    the   determination  of   the   outgoing  arc,    the   use   of   depth   in 
the   update   of   the   simplex  multipliers   is   more    critical   to   the   success   of   GNET   than   its 
use   to   find   the   outgoing  arc.      We   do  not  know   the   pivot   details   or   the   use  of   depth 
(if   any)    in   other   codes  nor  do  we  know   if   any  have   a  one-pass   update. 

Srinivasan   and  Thompson    [54]    1972  have  proposed  number  of  successors   to  be   used 
with   depth   to   determine  how  many  nodes   are   in   the   subarborescence   that   is    rehung  in 
the  pivot.      If   the   subarborescence  has   less    than  half  the  nodes   of   the   problem,    they 
propose  moving   the   root   and  performing   the  pivot   on    the   smaller  part  of   the  prede- 
cessor graph.      Our  experiments   show   that   the   subarborescence   can   always  be  expected   to 
have   significantly   less    than  half   the  nodes,    so  we  have   rejected   the   idea  of  moving 
the   root.      However,   we  have   discovered   that  number  of  successors   can  be   used  in  a  way 
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not  evisaged  by   Srinivasan  and  Thompson.      As    described  above   it   can   replace   depth   to 
support   the   determination   of   the   outgoing  arc  and  the  one-pass   update   of   the   simplex 
multipliers.      We  have   also  noted   that   the  number  of  successors   for  nodes   in   the  pivot 
subarborescence   changes   only   for  nodes   on   the   pivot  stem  avoiding  the   update   indicated 
in    [55].      Subsequent   to  our  presentation   of  GNET/Successors ,   Glover  and  Klingman    [28] 
19  75   have   proposed   a  number -of -successors   version   of  their  algorithm.      They   describe 
a  one-pass   update   of   the   simplex  multipliers    that     uses   an   additional  node   length 
array      to   store        for  each  node   the    last  successor   ranked  by  preorder.      They   conjec- 
ture   [28,   p.    7]    that   a   code  based  on   this   approach  will   dominate   their  previous    codes. 
This   may  be    true    for   their   codes,   but   as   discussed  above   this   is  not   our  experience 
in   comparing   our  one-pass    depth  version   against  our  one-pass  number  of  successors 
version.      Also, as   indicated  in   Figure    7, GNET/Successors   accomplishes   a  one-pass   update 
without   an   additional  array   to  store   and  manipulate. 

Preorder  distance  was  used  by  McBride  [44]  1973  and  Graves  and  McBride  [32]  1973. 
They  developed  the  zero- to-right  property  and  established  the  proposition  stated  above 
to    find    the   join. 

Mulvey    [45]    uses   a   candidate   list   of   arcs   that   corresponds    to   the   "partial 
suboptimization"    that   is   used   in   commercial   linear  programming  systems    [47].      His 
candidate    list   is    controlled  by   two  parameters:      x        and     x? .      In   our  terminology, 
general  scans   of  nodes   are    done   until    there   are      x„      arcs   on   the   candidate   list, 
then   up    to      x        pivots    are    performed   by    choosing   arcs    from   the    candidate    list.      The 
candidate    list    is    then   discarded  and   the   process  begins   again.      As  noted  in    [10] , 
Mulvey 's    approach   has   been    adopted  by   Glover   and  Klingman. 

The    candidate   queue    used   in   GNET   has   evolved    from  similar  mechanisms    developed 
by   Graves    for  more   general  mathematical   programming  systems.      It    is    unique    in    that   it 
contains  both  nodes   and  arcs,    it   is   used   cyclically   and   arcs   are  never   removed  as   long 
as    they    continue    to   price    out    favorably. 
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The  Srinivasan  and  Thompson  code  {55]  and  the  Harris  code  [33J  are  for  dense 

uncapacitated  transportation  problems;  all  other  contemporary  codes  known  to  the  authx 
solve  the  capacitated  transshipment  problem  with  a  sparse  representation  of  the  net- 
work similar  to  that  described  here. 

Harris  [33]  1976  has  described  a  code  for  dense  uncapacitated  transportation 

problems  that  have  few  sources  and  many  sinks.   The  brief  description  indicates  that 
the  simplex  multipliers  are  not  stored  for  the  sinks  but  there  are  no  details  on  the ■ 
handling  of  the  pivot  nor  on  the  use  of  a  traversal  mechanism.   Independently,  the 
authors  have  developed  TNET  for  the  sparse  capacitated  transportation  problem  and 
subsequently  developed  XNET  for  the  general  capacitated  transshipment  problem.   TNET 
and  XNET  are  specific  examples  of  the  generation  of  simplex  multipliers;  other 
refinements  that  generate  more  (or  all)  the  multipliers  are  possible. 

There  is  considerable  literature  on  postoptimal  analysis  and  parametric 
programming  for  transportation  and  transshipment  problems  (e.g.  [9,  52,  53]). 

CONCLUSION 

The  GNET  programs  are  small,  fast  and  easy  to  modify.   They  integrate  many  well- 
known  techniques  from  mathematical  optimization  and  computer  science.   It  is  importan 
however,  to  carefully  discriminate  between  the  underlying  ideas  of  mathematical  progr 
ming  and  the  computer  science  topics  applied  to  their  implementation.   A  sound 
mathematical  footing  is  required  for  generalization  beyond  pure  networks.   The 
computer  science  literature  has  contributed  a  great  deal  to  the  local  representation 
of  global  information  in  trees  and  graphs,  and  has  given  valuable  recursive  methods  f 
manipulating  data  structures;  many  of  these  techniques  can  be  applied  to  the  basic 
arborescence.   But  some,  such  as  rerooting,  pruning,  balancing  and  even  conversion  to 
an  equivalent  binary  tree  with  extra  dummy  nodes  and  arcs,  do  not  work  well  for  netwo 
problems.   In  these  cases,  the  mathematical  interpretation  implies  (as  does  the  exper 
mental  evidence)  that  these  devices  are  needless  complications  that  increase  solution 
expense  or  introduce  superfluous  equations  and  variables. 
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In  our  experience  large  scale  problems  are  always  created  from  source  data  by 
problem  generator  programs.   The  problem  generator  may  occupy  significantly  less 
computer  storage  than  the  file  of  coefficients  produced.   Thus,  it  can  be  worthwhile 
to  generate  the  data  as  needed  rather  than  store  it  explicitly.   Such  generators  are 
easily  incorporated  in  GNET.   (This  approach  does  not  entirely  avoid  arc  length  infor- 
mation for  capacitated  problems  since  a  record  must  be  maintained  for  each  arc  out  of 
the  basis  at  its  upper  bound.)   These  advantages  also  invite  development  and  use  of  a 
random  problem  generator  subroutine  to  replace  the  cumbersome  problem  files  and  cost 
of  using  NETGEN  [39]  to  generate  very  large  problems.   The  design  of  GNET  is  also 
consistent  with  the  use  of  explicit  arc  arrays  stored  in  peripheral  devices. 

Perhaps  the  most  important  potential  for  the  pure  network  solution  speed  of  GNET 
lies  in  more  general  large  scale  models  with  imbedded  networks  in  their  structure. 
The  multicommodity  distribution  system  design  code  of  Geoff rion  and  Graves  [22]  has 
been  revised  to  incorporate  GNET  to  repeatedly  solve  the  (thousands  of)  network  sub- 
problems.   A  goal  programming  code  for  network  problems  with  quadratic  objective 
functions  has  been  successfully  built  by  the  authors  with  GNET  used  as  the  key  sub- 
routine. 

The  theme  of  GNET  is  replacement  of  arithmetic  primal  simplex  computations  by 

simple  but  equivalent  logical  tests.   This  is  reminiscent  of  the  motives  for  generalized 
upper  bounding  and  suggests  that  network  factorization  may  prove  to  be  a  viable  com- 
petitor for  GUB  in  general  linear  programming  systems  (31,  32). 

Lee  [43]  indicates  that  truly  huge  models  may  be  solved  by  mathematical  aggrega- 
tion and  develops  a  wide  class  of  network  aggregation  methods  producing  surrogate 
problems  that  are  pure  networks  which  are  smaller  in  size  and  which  preserve  and 
exploit  special  global  structure  in  the  original  problem.   This  approach  is  intriguing 
due  to  the  curious  general  improvement  in  performance  encountered  when  solving  real 
models  rather  than  random  test  problems  of  equivalent  size. 
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The  FORTRAN  program  GNET/Depth  [6]    1975  is  distributed  for  a  nominal  handling 
charge  on  an  exclusive  use  basis.   For  further  information  write  Prof.  Glenn  Graves  a 
the  Western  Management  Science  Institute,  Graduate  School  of  Management,  University 
of  California,  Los  Angeles,  California  90024. 


54 


REFERENCES 

[1]      Balas,   E.    and  Hammer,   P.    L. ,    "On   the  Transportation  Problem  -  Part   I,"   Cahiers 
du  Centre   d' Etudes   de   Recherche   Operationelle ,    4^   No.    2,    1962,   p.    98. 

[2]      Barr,    R.    S.,    Glover,    F.    and  Klingman ,   D. ,    "An   Improved  Version  of  the  Out-of- 
kilter  Method  and  a  Comparative  Study   of   Computer  Codes,"  Mathematical 
Programming,  ]_,   No.    1,   August    1974,    p.    60. 

[3]      Bartels,    R.    H.    and  Golub ,    G.   H. ,    "The   Simplex  Method  of  Linear  Programming  Using 
LU  Decomposition,"   Communications   of   the  Association   for  Computing  Machinery,    12, 
No.    5,   May   1970,    p.    266. 

[4]      Berge,    C. ,    The  Theory   of   Graphs   and  its   Applications,    translated  by  A.    Doig, 
John  Wiley   and  Sons,   New  York,    1962. 

{5]      Bradley,    G.    H.  ,    "Survey   of  Deterministic  Networks,"  AIIE  Transactions,  _7,   No.    3, 
September   19  75,    p.    222. 

[6]      Bradley,    G.    H. ,   Brown,    G.    G.    and  Graves,    G.   W. ,    "GNET,   A  Primal  Capacitated 
Network  Program."      Copyright   1975. 

[7]      Busacker,    R.    G.    and  Gowen,    P.    J.,    "A  Procedure    for  Determining  a  Family   of 

Minimum-Cost  Network   Flow  Patterns,"   PRO  Technical  Report   15,    Operations   Research 
Office,   Johns   Hopkins   University,    1961. 

[8]      Busacker,    R.    G.    and  Saaty,   T. ,    Finite  Graphs   and  Networks ,   McGraw  Hill,   New  York, 
19  70. 

[9]      Charnes ,    A.    and  Cooper,   W.   W. ,   Management  Models   and  Industrial  Applications   of 
Linear  Programming,   Volumes   I   and  II,   John  Wiley   and  Sons,   New  York,    1961. 

[10]      Charnes,    A.,    Glover,    F. ,   Karney,    D. ,   Klingman,   D.    and  Stutz,   J.,    "Past,   Present 
and  Future   of   Large   Scale  Transshipment   Computer  Codes   and  Applications," 
Research   Report    CS    131,    Center   for   Cybernetic  Studies,   The   University   of  Texas, 
Austin,   July    1973    (revised  October   1974). 

[11]      Cheong,   Y.    P.,    "Network  Transformations   and  Some  Applications,"  MS  Thesis, 
Naval  Postgraduate   School,   December   1975. 

[12]      Clasen,    R.    J.,    "The  Numerical  Solution   of  Network  Problems   Using  the  Out-of-kilter 
Algorithm,"   RAND  Memorandum  RM-5456-PR,    Santa  Monica,    California,   March   1968. 

[13]      Dantzig,    G.    B. ,    "Application   of  the   Simplex  Method  to   a  Transportation  Problem," 
in  Activity  Analysis   of  Production   and  Allocation,   edited  by  T.    C.    Koopmans , 
John  Wiley   and  Sons,   New  York,    1951. 

[14]      Dantzig,    G.    B. ,    Linear  Programming  and  Extensions,   Princeton  University  Press, 
Princeton,   New  Jersey,    1963. 

[15]  Dantzig,  G.  B.  and  Van  Slyke,  R.  M. ,  "Generalized  Upper  Bounding  Techniques  for 
Linear  Programming,"  Journal  of  Computer  and  System  Sciences,  J.,  No.  3,  October 
1967,   p.    213. 


55 


Edmonds,   J.    and  Karp ,    R.    M. ,    "Theoretical  Improvements   in  Algorithmic  Efficiency 
for  Network  Flow  Problems,"  Journal  of  the  Association   for   Computing  Machinery, 
19,   No.    2,   April  1972,   p.    248. 

Elmaghraby,    S.,    Some  Network  Models   in   Operations   Research,    Springe r-Verlag, 
New  York,    19  70. 

Ford,    L.    R.    and  Fulkerson,    D.    R. ,    "A  Primal-Dual  Algorithm  for  the   Capacitated 
Hitchcock  Problem,"  Naval  Research  Logistics   Quarterly,    h_,  No.    1,  March   1957, 
p.    47. 

Ford,    L.    R.    and   Fulkerson,    D.    R. ,    Flows   in  Networks ,   Princeton  University  Press, 
Princeton,   New  Jersey,    1962. 

Fulkerson,    D.    R. ,    "An  Out-of-Kilter  Method  for  Minimal-Cost   Flow  Problems,"   SI AM 
Journal  of  Applied  Mathematics,    9_,   No.    1,   March   1961,   p.    18. 

Fulkerson,  D.  R. ,  "Flow  Networks  and  Combinatorial  Operations  Research,"  American 
Mathematical  Monthly,    73,   No.    2,    February   1966,   p.    115. 

Geoffrion,  A.  M.  and  Graves,  G.  W. ,  "Multicommodity  Distribution  System  Design  by 
Benders   Decomposition,"  Management   Science,    20,   No.    5,   January   1974,   p.    822. 

Gill,  P.  E.  and  Murray,  W. ,  "A  Numerically  Stable  Form  of  the  Simplex  Algorithm," 
Journal  of   Linear  Algebra  and  its   Applications,"   6_,    1973,   p.    99. 

Gilsinn,  J.  and  Witzgall, C. ,  "A  Performance  Comparison  of  Labeling  Algorithms  for 
Calculating  Shortest  Path  Trees,"  Technical  Note  772,  National  Bureau  of  Standard 
Washington,    D.    C.  ,   May   1973. 

Glicksman,    S.,   Johnson,    L.    and  Eselson,    L. ,    "Coding   the  Transportation  Problem," 
Naval  Research   Logistics   Quarterly,    ~1_,   No.    2,   June    1960,    p.    169. 

Glover,    F. ,   Karney ,    D.    and  Klingman,    D. ,    "Implementation   and  Computational 
Comparisons   of  Primal,    Dual   and  Primal-Dual   Computer   Codes    for  Minimum  Cost 
Network   Flow  Problems,"  Networks,    U,   No.    3,    19  74,   p.    191. 

Glover,    F. ,   Karney,   D. ,   Klingman,    D.    and  Napier,   A.,    "A  Computational  Study   on 
Start  Procedures,   Basis    Change   Criteria  and  Solution  Algorithms    for  Transporta- 
tion Problems,"  Management   Science,    20,   No.    5,   January    1974,   p.    793. 

Glover,    F.    and  Klingman,    D. ,    "Improved   Labeling   of   L.    P.    Bases    in  Networks," 
Research   Report   CCS   218,    Center   for   Cybernetic  Studies,   The   University   of  Texas, 
Austin,    August    1974    (revised  October   1975). 

Glover,    F. ,   Klingman,    D.    and  Stutz,   J.,    "Augmented  Threaded  Index  Method   for 
Network  Optimization,"   INFOR,    12,   No.    3,    October   1974,   p.    293. 

Graves,    G.    W. ,    "A  Complete   Constructive   Algorithm  for  the  General  Mixed  Linear 
Programming  Problem,"  Naval  Research   Logistics   Quarterly,    12 ,   No.    1,   March   1965, 
p.    1. 

[31]      Graves,    G.    W.    and  McBride,    R.    D. ,    "The   Factorization  Approach   to   Large-Scale 

Linear  Programming,"  Eighth   International  Symposium  on  Mathematical  Programming, 
Stanford  University,   August    1973    (also  in  Mathematical  Programming,    10,   No.    1, 
February   1976,   p.    91)  . 

56 


32]      Graves,   G.    W.    and  McBride,   R.    D.  ,    "A  Factorization  Algorithm  for  Network  Problems 
with  Side  Constraints,"   44th  National  ORS A  Meeting,    San  Diego,    Fall  1973. 

33]      Harris,   B. ,    "A  Code   for  the  Transportation  Problem  of  Linear  Programming," 

Journal  of   the  Association   for   Computing  Machinery,    23,   No.    1,   January   1976, 
p.    155. 

34]      Harris,    C.    C.    and  Hopkins,    F.    E. ,    Locational  Analysis   -  An  Interregional 

Econometric  Model   of  Agriculture,   Mining,  Manufacturing,  and  Services,    Lexington 
Books,    Lexington,   Massachusetts,   1972. 

35]      Hatch,    R.    S.,    "Bench  Marks   Comparing  Transportation   Codes  based  on  Primal  Simplex 
and  Primal-Dual   Algorithms,"   Operations   Research,    23,   No.    6,   November   1975, 
p.    1167. 

36]      Hitchcock,    F.    L. ,    "The   Distribution  of   a  Product    from  Several  Sources   to  Numerous 
Localities,"   Journal  of  Mathematics   and  Physics,    20,   No.    2,   April   1941,    p.    224. 

37]      Johnson,   E.    L. ,    "Networks   and  Basic  Solutions,"   Operations   Research,    14,   No.    4, 
August   1966,   p.    619. 

38]      Klein,   M. ,    "A  Primal  Method   for  Minimal   Cost   Flows  with  Application   to  the 
Assignment    and  Transportation  Problems,"  Management   Science,    14,   No.    3, 
November   196  7,   p.    205. 

39]      Klingman,    D. ,   Napier,    A.    and  Stutz,   J.,    "NETGEN  -  A  Program  for  Generating  Large 
Scale    (Un)    Capacitated  Assignment,   Transportation   and  Minimum  Cost   Flow  Network 
Problems,"   Management   Science,    20,   No.    5,   January   1974,   p.    814. 

40]      Koopmans ,   T.    C. ,    "Optimum  Utilization  of  the  Transportation  System,"  Proceedings 
of   the   International  Statistical  Conferences,   Washington,    D.    C. ,    1947,   published 
in  Volume   5,    1949,    p.    136.       (Also  in  Scientific  Papers   of  Tj ailing   C.    Koopmans, 
Springer- Verlag,   New  York,    19  70,   p.    184. 

41]      Knuth,    D.    E.,    The   Art    of   Computer  Programming,    Volume    1    (Fundamental  Algorithms), 
Addis on-Wes ley ,    Reading,   Massachusetts,    1968. 

42]      Langley,    R.    W. ,   Kennington,    J.    and  Shetty,    C.    M. ,    "Efficient    Computational 

Devices    for  the   Capacitated  Transportation  Problem,"  Naval  Research   Logistics 
Quarterly,    21,   No.    4,    December   1974,   p.    637. 

43]      Lee,    S.,    "Surrogate   Programming  by  Aggregation,"  Ph.D.    dissertation,   UCLA, 
September   19  75. 

44]      McBride,    R. ,    "Factorization   in   Large-Scale   Linear  Programming,"  Working  Paper 
No.    200    (also  Ph.D.    dissertation),   Western  Management   Science   Institute,   UCLA, 
June   19  73. 

45]      Mulvey,   J.,    "Column  Weighting   Factors   and  Other  Enhancements   to  the   Augmented 
Threaded  Index  Method   for  Network  Optimization,"  Joint  National  Meeting  of 
ORSA/TIMS,    San   Juan,    Puerto   Rico,    Fall   19  74. 

46]      Mulvey,   J.,    "Special  Structure   in  Network  Models   and  Associated  Applications," 
Ph.D.    dissertation,   UCLA,   August    19  75. 


57 


[47]      Orchard-Hays,   W. ,   Advanced  Linear-Programming   Computing  Techniques ,   McGraw-Hill, 
New  York,    1968. 

[48]      Orden,   A.,    "The  Transshipment   Problem,"   Management   Science,   2,   No.    3,   April   1956. 

[49]      Saunders,   M.    A.,    "Large   Scale   Linear  Programming  using  the   Cholesky   Factorization 
Technical  Report   CS-72-252,    Computer  Science   Department,    Stanford  University,    197 

[50]      Saunders,   M.    A.,    "Product   Form  of   the   Cholesky   Factorization   for  Large-Scale 
Linear  Programming,"  Technical  Report   CS-72-301,    Computer  Science   Department, 
Stanford  University,    1972. 

[51]      Scoins,   H.    I.,    "The    Compact   Representation   of  a  Rooted  Tree   and  the  Transportatic 
Problem,"   International  Symposium  on  Mathematical  Programming,    London,    1964. 

[52]      Srinivasan,   V.    and  Thompson,    G.    L. ,    "An  Operator  Theory   for  the  Transportation 
Problem  -   I,"  Naval  Research   Logistics   Quarterly,    19 ,   No.    2,   June   1972,   p.    205. 

[53]      Srinivasan,    V.    and  Thompson,   G.    L. ,    "An   Operator  Theory   for  the  Transportation 
Problem  -   II,"  Naval  Research   Logistics   Quarterly,    19,   No.    2,   June   1972,   p.    227. 

[54]      Srinivasan,   V.    and  Thompson,    G.    L. ,    "Accelerated  Algorithms   for  Labeling   and 
Relabeling   of  Trees  with  Application   for  Distribution  Problems,"  Journal  of 
the   Association   for  Computing  Machinery,    19 ,   No.    4,    October   1972,   p.    712. 

[55]      Srinivasan,   V.    and  Thompson,    G.    L. ,    "Benefit-Cost   Analysis   of  Coding  Techniques 
for  the  Primal  Transportation  Algorithm,"  Journal   of  the  Association   for 
Computing  Machinery,    20,   No.    2,    April   1973,   p.    194. 


58 


INITIAL  DISTRIBUTION  LIST 


Copies 


Defense  Documentation  Center  2 

Cameron  Station 
Alexandria,  Virginia   22314 

Library,  Code  0212  2 

Naval  Postgraduate  School 
Monterey,  California  93940 

Dean  of  Research,  Code  023  2 

Naval  Postgraduate  School 
Monterey,  California  93940 

Library,  Code  55  1 

Naval  Postgraduate  School 
Monterey,  California  93940 


59 


11 1<  o?5l 


DUDLEY  KNOX  LIBRARY  -  RESEARCH  REPORTS 


5  6853  01071684  8 


U176 


f\/P^ 


